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Abstract 

This article reviews recent developments in the theoretical under- 
standing and the numerical implementation of variational renormalization 
group methods using matrix product states and projected entangled pair 
states. 
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One of the biggest challenges in physics is to develop efficient methods for 
simulating Hamiltonians: the basic laws of quantum mechanics are very well 
understood, but the crucial understanding of the collective behavior of inter- 
acting elementary particles and of emerging phenomena relies on our ability to 
approximate or simulate the corresponding equations. Simulation is particularly 
relevant in the context of strongly interacting quantum systems, where conven- 
tional perturbation theory fails. Popular methods to simulate such systems 
are the numerical renormalization group algorithms by K. Wilson |140| and S. 
White [T33H31GI71G1I] and quantum Monte Carlo methods. Both of these have 
been used with great success, but they also have severe limitations: the numer- 
ical renormalization group and DMRG only work in cases where the original 
Hamiltonian can be mapped to a local Hamiltonian defined on a 1-dimensional 
chain, while Monte Carlo methods suffer from the so-called sign problem [113j 
which makes them inappropriate for the description of fermionic and frustrated 
quantum systems. Very recently however, new insights coming from the field of 
quantum information and entanglement theory have shown how those numerical 
RG-methods can be generalized to higher dimensions, to finite temperature, to 
random systems etc. [UH HH EH EZU EH EH E2 EH M E] These develop- 
ments bear the promise to revolutionize the way quantum many-body systems 
can be simulated and understood, which will prove essential in the light of the 
ongoing miniaturization of electronic devices and of fundamental problems such 
as the study of the role of the Hubbard model in high T c superconductivity. 

In this review, we will focus on those new developments. In particular, we 
will focus on the simulation of quantum lattice systems of spins and/or fermions. 
The Hamiltonians arising in this context can be thought of as effective Hamil- 
tonians of more complicated systems that capture the physics for the relevant 
low-energy degrees of freedom. Studying quantum spin systems has a long and 
very interesting history, starting with Dirac and Heisenberg in the 1920's where 
they proposed the so-called Heisenberg model [SUJ HHj as being illustrative for 
the basic mechanism giving rise to magnetism. The simulation of quantum spin 
systems and fermionic systems on a lattice has however turned out to be ex- 
tremely complicated. Consider for example the Hubbard model in 2 dimensions 
that is representing the simplest possible model of fermions hopping on a lat- 
tice and exhibiting on-site interactions: despite considerable efforts because of 
the connection with high T c superconductivity, it is absolutely unclear how the 
low-energy wavefunctions look like in the relevant parameter range. This is the 
most notable illustration of the biggest bottleneck in condensed matter theory: 
numerical simulation of effective simple Hamiltonians. The new tools described 
in this review bear the promise of opening up that bottleneck. 

This review is organized as follows. In a first, rather technical, section, we 
will review basic properties of quantum spin systems and provide the intuitive 
justification for the numerical tools used later. This section can be skipped 
by those who are more interested in the variational numerical renormalization 
methods and less in the underlying theory; however one of the major advantages 
of the methods is precisely that there is a solid theoretical justification for using 
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them. In section [3J we review the numerical renormalization group method as 
introduced by Wilson, and show how this gives rise to the concept of matrix 
product states. Section 3 discusses the features of matrix product states, and 
goes on to applications on spin chains such as a reformulation of the density 
matrix renormalization group (DMRG) and generalizations that allow to treat 
excited states. Section 2] is devoted to the issue of real and imaginary time- 
evolution of 1-dimensional quantum spin systems, and section[5]to the concept of 
matrix product operators leading to extensions of DMRG to finite temperature 
systems, to random spin chains and to the simulation of 2-D classical partition 
functions. In section we will show how all these ideas can be generalized 
to the 2-dimensional case by introducing the class of Projected Entangled Pair 
States (PEPS), and we will discuss the applicability and limitations of these 
methods. Finally, we discuss the convex set of reduced density operators of 
translational invariant states in appendix El the relation between block entropy 
and the accuracy of a matrix product state approximation in appendix [Bj and 
we include some explicit matlab programs of matrix product state algorithms 
in appendix [C] 

At this point, we would like to warn the reader that this is not a review 
of DMRG methods and possible extensions. In fact, there exist very good 
reviews of that topic (see, for example, [97l [90] ) , where the important progress 
experienced by DMRG-like methods is extensively discussed. Here, we offer 
a review of the new methods that have been introduced during the last few 
years, and that have evolved from ideas in the context of Quantum Information 
Theory. To make the review self-contained and uniform, we will mostly focus on 
the methods introduced by ourselves and collaborators, and will also touch upon 
the approach advocated by Vidal. The main reason is that they can be viewed 
as (generalized) variational techniques over sets of entangled states, and thus we 
can present in a unified form the vast variety of algorithms whose description 
is thusfar scattered in several publications. Furthermore, we will also discuss 
some issues which have not be published so far. 
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1 Spin systems: general features 



One of the main characteristics of quantum mechanics is that the underlying 
Hilbert space is endowed with a tensor product structure: the Hilbcrt space 
of two interacting systems is given by the tensor product space of the two 
individual ones. This structure of the Hilbert space is a direct consequence 
of the superposition principle, and opens up the possibility of entanglement 
and new phases of matter. This simple tensor product structure, however, 
makes it very clear what the main obstacle is in developing a general theory of 
many-body quantum systems: the size of the Hilbert space grows exponentially 
in the number of basis constituents, and hence we would, in principle, need 
exponentially many variables to specify the wavefunction of a iV-particle system. 
However, it is a priori not clear whether Nature can fully exploit and explore 
these vast territories of Hilbert space, because another main characteristic of 
quantum mechanics is that interactions always seem to happen locally and only 
between a few bodies. 

As it turns out, all physical states live on a tiny submanifold of Hilbcrt 
space 0- This opens up a very interesting perspective in the context of the 
description of quantum many-body systems, as there might exist an efficient 
parametrization of such a submanifold that would provide the natural language 
for describing those systems. In the context of many-body quantum physics, 
one is furthermore mainly interested in describing the low energy sector of local 
Hamiltonians, and as we will discuss later, this puts many extra constraints on 
the allowed wavef unctions. 

As a trivial example of such a submanifold, let us consider a system of 
N — > oo spins that all interact with each other via a permutation invariant 2- 
body Hamiltonian such as the Heisenberg interaction J5UJ US] ■ Note that there 
is always a ground state that exhibits the same symmetry as the Hamiltonian. 
As a consequence of the quantum de-Finnetti theorem 108J one can show that 
the manifold of all density operators with permutation symmetry is exactly 
the set of all separable states p = ^2iPipf N when N — > oo. Ground states 
correspond to the extreme points of this set, which are exactly the product 
states that have no entanglement. All ground states of permutational invariant 
systems therefore lie on the submanifold of separable states which indeed have 

J Let us for example consider a quantum system of N spin 1/2's on a cubic lattice and 
all spins polarized in the z-direction, and let us pick a random state in the 2^ dimensional 
Hilbert space according to the unitarily invariant Haar measure. Let us furthermore ask the 
question to find a lower bound on the time it would take to evolve the polarized state into one 
that has an overlap with the random one that is not exponentially small, and this with local 
interactions. Using tools developed in the context of quantum information theory, one can 
show that this lower bound scales exponentially in the number of spins [123] : any evolution 
over a time that scales only polynomially in the number of spins would only allow to get an 
exponentially small overlap with the chosen random state. If we would for example consider 
a system of a few hundred of spins, it would already take much longer than the lifetime of 
the universe to come close to any random point in its Hilbert space. This shows that almost 
all points in the Hilbert space of a many-body quantum system are unphysical as they are 
inaccessible: all physical states, i.e. states that can ever be created, live on a tiny submanifold 
of measure zero. 
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an efficient representation (only N vectors, one for each individual spin, have to 
be specified); this is the equivalent statement as saying that mean- field theory 
becomes exact in the thermodynamic limit. This can also be understood in the 
context of the monogamy property of entanglement [TH [S3] : a spin has only a 
finite entanglement susceptibility, and if it has to share the entanglement with 
infinitely many other spins then the amount of entanglement between 2 spins 
will be go to zcrod. 

In the more general situation when no permutation symmetry is present 
but only a smaller symmetry group such as the group of translations, the char- 
acterization of the relevant manifold is much harder. In contrast to the case 
of permutation symmetry, where every pure state of N spin s systems can be 
written as a linear combination of at most N 2s /(2s)\ Dicke states 0, the size of 
a complete basis of translational invariant states is exponentially increasing as 
(2s + l) N /N, and hence the manifold of interest has exponentially many pa- 
rameters. But ground states of local Hamiltonians have many more nongencric 
properties, most notably the fact that they have extremal local properties such 
as to minimize the energy: as the energy is only dependent on the local prop- 
erties, and the ground state is determined by the condition that its energy is 
extremal, ground states have extremal local properties and the global properties 
only emerge to allow for these local properties to be extremal. As an example, 
let us consider a spin 1/2 antiferromagnetic chain with associated Hamiltonian 



where the notation denotes the sum over nearest neighbours. Each in- 

dividual term in the Hamiltonian corresponds to an exchange interaction and 
would be minimized if spins i and j are in the singlet state 



but due to the monogamy or frustration properties of entanglement, a spin 1/2 
cannot be in a singlet state with more than one neighbour. As a result, the 
ground state becomes a complicated superposition of all possible singlet cover- 
ings, and as an interesting by-product quasi long-range order may arise. The 
important point, however, is that this wavefunction arises from the condition 
that its local properties are extremal: finding ground states of local translational 
invariant 2-body Hamiltonians is equivalent to characterizing the convex set of 

2 In the case of a translational and rotational invariant system on a regular lattice, the finite 
versions of the de-Finnctti theorem 61, 16 allow one to find lower bounds on the distance 
of the reduced 2-body density operators to the separable ones; a scaling of the form 1/c is 
obtained where c is the coordination number. See also appendix A. 

3 Dicke states are well studied in the context of quantum optics, and are obtained by taking 
the linear superposition of all permutations of a pure separable state in the standard basis; in 
the case of spin 1/2, you can label them by their expectation value of Sf. See e.g. 11071 
for an extensive discussion of their entanglement properties. 
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2-body density operators compatible with the fact that they originate from a 
state with the right symmetry (see appendix [X]) . 

Clearly, we would like to parameterize the manifold of states {iV'ea}} with 
extremal local properties. In practice, it is enough to parameterize a manifold 
of states {\ipappr}} such that there always exist a large overlap with the exact 
states {\ipex}}- V\">Pex)i3\ipappr} ■ \\\i>ex) ~ \^ap P r)\\ < £• Let us consider any 
local Hamiltonian of N spins that exhibits the property that there is a unique 
ground state \tpex) an d that the gap is A (AT). Let us furthermore consider the 
case when A (AT) decays not faster than an inverse polynomial in N (this con- 
dition is satisfied for all gapped systems and for all known critical translational 
invariant systems in ID). Then let us assume that there exists a state \ipappr) 
that reproduces well the local properties of all nearest neighbour reduced den- 
sity operators: || p ap pr — Pex\\ < S. Then it follows that the global overlap is 
bounded by 

HIVw) - IVw)ll 2 < AfM' 

This is remarkable as it shows that it is enough to reproduce the local properties 
well to guarantee that also the global properties are reproduced accurately: for a 
constant global accuracy e, it is enough to reproduce the local properties well to 
an accuracy S that scales as an inverse polynomial (as opposed to exponential) in 
the number of spins. This is very relevant in the context of variational simulation 
methods: if the energy is well reproduced and if the computational effort to get 
a better accuracy in the energy only scales polynomially in the number of spins, 
then a scalable numerical method can be constructed that reproduces all global 
properties well (here scalable means essentially a polynomial method). 

The central question is thus: is it possible to find an efficient parameteriza- 
tion of a manifold of states whose local properties approximate well all possible 
local properties? A very interesting new development, mainly originating in the 
field of quantum information and entanglement theory, has shown that this is 
indeed possible. The main idea is to come up with a class of variational wave 
functions that captures the physics of the low-energy sector of local quantum 
spin Hamiltonians. 

So what other properties do ground states of local quantum Hamiltonians 
exhibit besides the fact that all global properties follow from their local ones? 
The key concept to understand their structure is to look at the amount of en- 
tanglement present in those states [13 1J : entanglement is the crucial ingredient 
that forces quantum systems to behave differently than classical ones, and it is 
precisely the existence of entanglement that is responsible for such exotic phe- 
nomena like quantum phase transitions and topological quantum order [133U69] . 
It is also the central resource that gives rise to the power of quantum comput- 
ing [76], and it is known that a lot of entanglement is needed between the 
different qubits as otherwise the quantum computation can be simulated on a 
classical computer [57l 1126] . This is because the amount of entanglement ef- 
fectively quantifies the relevant number of degrees of freedom that have to be 
taken into account, and if this is small then the quantum computation could 
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Figure 1 : The entropy of a block of spins scales like the perimeter of the block. 

be efficiently simulated on a classical computer. In the case of ground states 
of strongly correlated quantum many-body systems, there is also lot's of entan- 
glement (in the case of pure states, the connected correlation functions can be 
nonzero iff there is entanglement), but the key question is obviously to ask how 
much entanglement is present there: maybe the amount of entanglement is not 
too big such that those systems can still be simulated classically? 

Let us for example consider a quantum spin system on a n-dimensional 
infinite lattice, and look at the reduced density operator p^ of a block of spins 
in an L x L x ...L hypercube (Figure [T]). The von-Neumann entropy of pl is a 
coarse-grained measure that quantifies the number of modes in that block that 
are entangled with the outside [76], and the relevant quantity is to study how 
this entropy scales with the size of the cube. This question was first studied in 
the context of black-hole entropy [SJ 106, 38, HI [35] and has recently attracted 
a lot of attention |131| [9TJ [T2] . Ground states of local Hamiltonians of spins or 
bosons seem to have the property that the entropy is not an extensive property 
but that the leading term in the entropy only scales as the boundary of the 
block (hence the name area- law): 

S( PL ) ~ cL n -\ (1) 

This has a very interesting physical meaning: it shows that most of the entan- 
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Figure 2: A one dimensional spin chain with finite correlation length £ corr ; Iab 
denotes the distance between the block A (left) and B (right). Because Iab 
is much larger than the correlation length £ corr , the state pab is essentially a 
product state. 

glement must be concentrated around the boundary, and therefore there is much 
less entanglement than would be present in a random quantum state (where the 
entropy would be extensive and scale like L n ). This is very encouraging, as it 
indicates that the wavcfunctions involved exhibit some form of locality, and we 
might be able to exploit that to come up with efficient local parametcrizations 
of those ground states. 

The area law ([1]) is mildly violated in the case of 1-D critical spin systems 
where the entropy of a block of spins scales like [1311 [T2"] 

<?( Pi )^£±£io g L, 

but even in that case the amount of entanglement is still exponentially smaller 
then the amount present in a random state. It is at present not clear to what 
extent such a logarithmic correction will occur in the case of higher dimensional 
systems: the block-entropy of a critical 2-D system of free fermions scales like 
L log L [1421 [5^1 US] , while critical 2-D spin systems were reported where no 
such logarithmic correction are present [124j . but in any case the amount of 
entanglement will be much smaller than for a random state. It is interesting 
to note that this violation of an area law is a pure quantum phenomenon as it 
occurs solely at zero temperature: in a recent paper [141] . it has been shown 
that the block entropy, as measured by the mutual information (which is the 
natural measure of correlations for mixed states), obeys an exact area law for all 
local classical and quantum Hamiltonians. The logarithmic corrections therefore 
solely arise due to the zero-temperature quantum fluctuations. From the prac- 
tical point of view, that might indicate that thermal states at low temperature 
are simpler to simulate than exact ground states. 

The existence of an area law for the scaling of entropy is intimately connected 
to the fact that typical quantum spin systems exhibit a finite correlation length. 
In fact, M. Hastings has recently proven that all connected correlation functions 
between two blocks in a gapped system have to decay exponentially as a function 
of the distance of the blocks [15]. Let us therefore consider a 1-D gapped 
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quantum spin system with correlation length ^ corr . Due to the finite correlation 
length, the reduced density operator pab obtained when tracing out a block C 
of length Iab 3> £corr (sec figure will be equal to 

Pab ^ Pa® Pb (2) 

up to exponentially small corrections 0- The original ground state \ipABc) is a 
purification of this mixed state, but it is of course also possible to find another 
purification of the form \ipACi) ® IV'scv) (up to exponentially small corrections) 
with no correlations whatsoever between A and B; here Ci and C r together span 
the original block C. It is however well known that all possible purifications of 
a mixed state are equivalent to each other up to local unitaries on the ancillary 
Hilbert space. This automatically implies that there exists a unitary operation 
Uc on the block C (see figure ^ that completely disentangles the left from the 
right part: 

1a ® U c <E> Ib\^abc) - \ipAd) ® IVtecv). 
This implies that there exists a tensor A l a * with indices 1 < a, (3, i < D (where 
D is the dimension of the Hilbert space of C) and states IV'a), iV'f ), defined 
on the Hilbert spaces belonging to A, B, C such that 

\i>ABG) * £ &MM)\tf)- 

Applying this argument recursively leads to a matrix product state (MPS) de- 
scription of the state (we will define those MPS later) and gives a strong hint 
that ground states of gapped Hamiltonians are well represented by MPS. It 
turns out that this is even true for critical systems [116j ; a proof is presented in 
appendix B. 

A remarkable feature of any gapped spin chain is thus that one could imagine 
dividing the whole chain in segments of size / ^> £ corT . , and then apply a disen- 
tangling operation on all blocks in parallel. This would then lead to a product 
state of many parts, and as such gives a procedure of how such a ground state 
could be prepared using a quantum circuit with a logical depth that is only 
dependent on the correlation length (and independent of the number of spins!). 

In the case of critical systems, we do not expect that this disentangling pro- 
cedure works as correlations on all length scales appear. However, G. Vidal 
showed how this can be remedied by introducing some more advanced disen- 
tangling scheme that acts on many different length scales [130] . Basically, the 
idea is follow up a disentangling step by a coarse-graining step, and do this 

4 Prom a purely mathematical point of view, this is not correct as, surprisingly, there exist 
states for which all connected correlations functions are negligable while they are very far from 
being tensor product states 1491 1481 : examples exist with negligable correlation functions but 
with the mutual information in the order of the number of qubits in every block. To remedy 
this, we pointed out in a recent paper that a more sensible way of defining a correlation 
length is by using the concept of mutual information, and when there is an exponential decay 
of mutual information, then the state pab is guaranteed to be close to a tensor product state 

OH- 
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recursively until there is only one spin left. The procedure for doing so is called 
the Multiscale Entanglement Renormalization Ansatz (MERA), and may lead 
to alternative methods for simulating quantum spin systems. 

It is interesting to note that the situation can again be very different in two 
dimensions: in that case, gapped quantum spin systems can exhibit topological 
quantum order, and in [5] it was proven that the depth of any quantum circuit 
preparing such a topological state has to scale linearly in the size of the system. 
In other words, it is impossible to device a scheme/pattern by which one could 
disentangle such states in a parallel way as can be done in 1-D. However, also 
in this 2-D case it is possible to come up with a variational class of states, the 
so-called projected entangled pair states |118l I75j . that captures the essential 
physics for describing those systems. By a recent argument of M. Hastings 
[171 H5] , it can again be proven that basically every ground state of a local 
Hamiltonian can be well represented by a state within this class. 
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2 Wilson's numerical renormalization group method 



The exact description of strongly correlated condensed matter systems poses 
formidable difficulties due to the exponentially large dimension of the associated 
Hilbert space. However, K. Wilson was the first one to understand that the 
locality of the interactions between particles or modes enforces ground states to 
be of a very specific form and could be exploited to simulate them. This insight 
led to the development of numerical renormalization group algorithms (NRG) 
|140| . The reason for the remarkable accuracy of NRG if applied to quantum 
impurity problems (e.g. the Kondo and Anderson Hamiltonians [65 , 166] ) can be 
traced back to the ability of mapping the related Hamiltonians to momentum 
space such that they become inhomogeneous 1-D quantum lattice Hamiltonians 
with nearest neighbor interactional. NRG is then a recursive method for finding 
the low-energy spectrum of such Hamiltonians, and yields very accurate results 
when there is a clear separation of energies, reflected by e.g. an exponential 
decay of the couplings within the 1-D hopping Hamiltonian. The NRG method 
then recursively diagonalizes the Hamiltonian from large to small energies: at 
each iteration, a tensor product of the larger energy modes with lower energy 
modes is made and then projected on a subspace of the lower energy modes 
of the combined system. Thereafter the Hamiltonian is rescaled. Hence the 
basic assumption is that the low energy modes are affected by their high energy 
counterparts, but not vice- versa. 

Let us illustrate Wilson's method on the hand of the NRG treatment of the 
impurity Anderson model (SIAM). This model can be mapped to a hopping 
Hamiltonian after a logarithmic discretization of the conduction band [66j : 



1 / U 



U = l^J n UV(«+Dm + h - c -) + D{ ed+ 2) C ^ Cdtl 

;(flc dfl + h.c.) + ^ (4 M c d „ - I)' (3) 

r-/ 2 q + r-^i-r-"- 1 ) 

2 v /( 1 _ r -2„-l)( 1 _ r -2„-3) 

Here /x can take the values J,, f, denotes the annihilation operator of the 
impurity and f nfjt of the n'th fermion with spin fi, summation over /i has been 
assumed, and N — * oo. As the dimension of the associated Hilbert space is 
2 2N , an exact diagonalization is impossible, and approximations must be made. 
The hopping terms are decaying exponentially in n, and the basic idea of NRG 
is to treat the largest no terms first which involves the diagonalization of a 
2 2n ° matrix. Next we specify a control parameter D < 2 2 ™ , retain only the 
eigenvectors 

{|C°>}a=l..D 



5 The basic reason why a 1-D model is obtained even though the original model is concerned 
with a magnetic impurity in a 3-D system is that the Kondo Hamiltonian only affects the s- 
wave part of the wavefunction; the s-wave modes are therefore dominant at low energies. 
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corresponding to the D lowest eigenvalues, and project the first no terms of the 
Hamiltonian onto that subspace using the projector 

p [n ° ] =X>xc°i 

a 

yielding the D x D matrix In the first step of the iteration, the extra 

term involving the coupling between the no'th and the (no + l)'th fermions is 
considered, and an exact diagonalization in the corresponding 4_D-dimcnsional 
Hilbert space is performed yielding the eigenvectors. To avoid that the dimen- 
sions of the effective Hamiltonian blow up, we project the Hamiltonian onto the 
new D-dimcnsional eigenspace corresponding to the lowest eigenvalues 

p["° +1 ] = £|«)(C 0+1 l 

a 

yielding 7i. no+1 . Note that p[™ + 1 l is a D x AD matrix, and for later reference 
we write its coefficients in tensor form as P l } n ^ +1 \ 1 < a, /3 < D, 1 < i < 4. 
Now we iterate this procedure N ~ n times, and NRG is typically said to have 
converged when T'H N ~ 1 = H N + est up to a unitary transformation. Note that 
the computational complexity of the NRG procedure scales as ND 3 . 

Let us now consider more closely the subspace on which we projected the 
original Hamiltonian. The D states {\ipa N }} & t the end of the iterations can be 
written as 

\lb N )~ V V pW"o]pin +i[no+l] piN [N] ,. v,. \ I- \ 

a no ...a nN _ 1 i„ in +l—i N 

(5) 

These states are exactly of the form as the ones mentioned in the introduction, 
and due to the feature that they are defined as a product of matrices, these 
are called matrix product states (MPS) [551(96, 88, 60, [HHj and were originally 
introduced in the mathematical physics community under the name of finitely 
correlated states [331 131] (a precursor of it appeared in the context of quantum 
Markov chains |Tj). The energies calculated using NRG are thus energies of 
an effective Hamiltonian which is the original one projected onto a subspace of 
MPS. In practice, the NRG method is highly successful for problems where the 
different terms in the Hamiltonian act on a different scale of energy, and in that 
case results up to essentially machine precision can be obtained; this can only 
be true if the class of the matrix product states indeed capture all the physics 
needed to describe the low-energy physics of these Hamiltonians. 

However, by looking at NRG by means of matrix product states, it is already 
clear that it can in principle be formulated as a variational method within the 
set of MPS. It turns out that this is exactly what S. White did by introducing 
the density matrix renormalization group (DMRG) [137| 1136] , but the way of 
looking at both methods from that point of view of MPS was only discovered 
much later [122j . 



13 



Singlet 



/ 

Projector 



Figure 3: Building up the AKLT state by partial projections on bipartite singlets 



3 Matrix product states and ground states of 
spin chains 

3.1 Construction and calculus of MPS 
3.1.1 The AKLT-model 

The notion of matrix product states (MPS) already appeared naturally in the 
section describing spin chains with finite correlation length and in the context of 
NRG. They were first studied in the work of Affleck. Kennedy, Lieb and Tasaki 
(AKLT) [2], where it was proven that the exact ground state of the spin-1 spin 
chain with Hamiltonian 

7~Laklt = ^ ySi.Sj + — ^Si.Sjj 



(id) 



3 



can be parameterized exactly as a matrix product state. To see this, they ob- 
served that the terms Pij are projectors {Pij) 2 = P%j onto the 5-dimensional spin 
2 subspace of 2 spin l's, and proceeded by constructing the unique ground state 
I^aklt) which is annihilated by all projectors Pij acting on nearest neighbours. 
This state \iPaklt) can be constructed as follows: 

• Imagine that the 3-dimcnsional Hilbcrt space of the a spin 1 particle is 
effectively the low-energy subspace of the Hilbcrt space spanned by 2 spin 
1/2's, i.e. the 3-D Hilbert space is the symmetric subspace of 2 spin 1/2 
particles. 

• To assure that the global state defined on the spin chain has spin zero, let 
us imagine that each one of the spin 1/2's is in a singlet state with a spin 
1/2 of its neighbours (see figure [3]). 

• The AKLT state can now be represented by locally projecting the pair of 
spin 1/2's in the symmetric subspace onto the spin-1 basis {|1), |0), | — 1)}: 

P = | i) AOOl-W i | 0) A Q1 l + ( 1Q h , hxAOOI + ^I 
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where t x = <j x , r y = ia y , t z = a z with {a a } the Pauli matrices and where 
we identified | - 1) = \x), |0) = \z), |1) = \y). 

Historically, the AKLT state was very important as it shed new lights into the 
conjecture due to Haldane [43] [42] that integer spin Heisenberg chains give rise 
to a gap in the spectrum. That is a feature shared by all generic matrix product 
states: they are always ground states of local gapped quantum Hamiltonians. 

Let us first try to rewrite \iPaklt) in the MPS representation. Let us assume 
that we have an AKLT system of N spins with periodic boundary conditions; 
projecting the wavefunction in the computational basis leads to the following 
identity: 

(0:1,0:2, —Oln\^AKLt) = Tr(r Ql .r y .r Q2 .T y ...r QN r y .) . 

The different weights can therefore be calculated as a trace of a product of 
matrices. The complete AKLT state can therefore be represented as 

\iPaklt) = ^2 Tr {T ai .Ty.T a2 .T y ...T a - N T y .) |oi)|q:2>---|on> 

which is almost exactly of the same form as the matrix product states introduced 
in [5] The only difference between them is the occurence of the matrices t v 
between the different products. This is however only a consequence of the fact 
that we connected the different nodes with singlets, and we could as well have 
used maximally entangled states of the form 

2 

|i)=£l«> 

i=l 

and absorbing r y into the projector; this way we recover the standard notation 
for MPS. 

So what did we learn from this AKLT-example? Basically, that there is a 
way of parameterizing the exact ground state of a particular strongly correlated 
quantum spin chain using a construction involving virtual bipartite entangle- 
ment and projections. From the point of view of quantum information theory, 
this parametrization is very appealing, as it gives an explicit way of construct- 
ing a highly entangled multipartite quantum state out of bipartite building 
blocks. More importantly, it is immediately clear how this picture can be gen- 
eralized [1201 1121] : instead of taking spin D = 2 bonds corresponding to spin 
1/2's, we can take much larger D for the virtual spins, and furthermore it is 
obvious that the projectors can be replaced with any linear map. What is really 
exciting in doing so is the fact that the states arising from this are transla- 
tional invariant by construction [33| : this is highly relevant as there does not 
seem to be another simple way of parameterizing translational invariant states. 
Recalling the discussion in section 1 , it is now obvious that this class of transla- 
tional invariant MPS , originally introduced in the literature under the name of 
finitely correlated states [33] , would form a very good ansatz for parameterizing 
ground states. Indeed, ground states of local Hamiltonians are characterized by 
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Map^i^i;i;s4«l B K*K'l 

n-l t-i 

Figure 4: A general Matrix Product State 

extremal local properties that are still compatible with the global translational 
symmetry, and as the MPS are build up by projecting the underlying maximally 
entangled states with extremal local correlations, it is very plausible that the 
MPS are perfectly suited for thifl 

Note that the number of parameters that we have to our disposition in these 
translational invariant MPS scales as dD 2 with d the physical dimension of the 
spins (indeed, we have to parameterize d matrices A a of dimension D x D). 
and hence the natural question to be asked now is how the convex set of all 
local reduced density operators so obtained compares to the exact convex set 
of all possible translational invariant systems. This problem is considered in 
Appendix [5] and [B] and it is found that the number of parameters needed has 
to scale as a constant or at most polynomially in the number of spins. 



3.1.2 Matrix Product States 

As already explained in the previous section, the obvious generalization of the 
AKLT-states is obtained by making the dimensions of the virtual spins D larger 
and to consider general linear maps A a instead of projectors, but we can make 
a further generalization by making D (the dimension of the virtual subsystems) 
and those maps site-dependent and write them as Di and A\ . The most general 
form of a MPS on N spins of dimension d is then given by 

IV') = E Tr(«..0|a 1 )|a 2 )...M (6) 

ai,Q2,---CKjV 

The matrices A % have dimension Di x Di+i (here we take the convention that 
Dn+i = Z?i), and a system with open boundary conditions is obtained by choos- 
ing D\ = 1 (i.e. no "singlet" between the endpoints). A pictorial representation 
is given in figure [4j Before continuing, let us remark that every state of N 
spins has an exact representation as a MPS if we let D to grow exponentially 
in the number of spins; this can easily be shown by making use of the tool of 

e Besides the AKLT- model, many variations of that Hamiltonian have been studied with ex- 
act MPS-states as ground states 64 62 , 63 ; this is particularly interesting because analytical 
solutions of spin chains are very rare. 
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Figure 5: a) An MPS with open boundary conditions as a tensor network. 
The open bonds correspond to the uncontractcd physical indices and the closed 
bonds to contracted indices arising from taking the products; b) calculating 
the expectation value of an nearest neighbour operator over a MPS with open 
boundary conditions by contracting the tensor network completely; c) an extra 
bond has to be contracted in the case of periodic boundary conditions. 



quantum teleportation as shown in |121j . However, the whole point of MPS is 
that ground states can typically be represented by MPS where the dimension D 
is small and scales at most polynomially in the number of spins; this is the ba- 
sic reason why renormalization group methods are exponentially more efficient 
than exact diagonalization. 



3.1.3 Calculus of MPS 

Let us next explain the calculus of those MPS. A crucial aspect of those MPS 
is indeed the fact that the expectation value of a large class of observables can 
be calculated efficiently. More specifically, this holds for all observables that 
are a tensor product of local observables, such as a l x <g> o{.. The easiest way to 
understand this is by considering a so-called tensor network, in which connected 
bonds correspond to contracted indices and open bonds to uncontracted ones 
[see figure [5] . More specifically, the tensor network in figure corresponds to 
the quantity 

A 1 A 2 A N 

as the "virtual" indices are all contracted (corresponding to taking products of 
matrices) and the physical ones are left open. To calculate expectation values 
of observables, it is clear that we also have to contract those physical indices, 
after sandwiching in the observables, and hence expectation values can be rep- 
resented by a completely contracted network. The relevant thing of course is 
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now to quantify the computational complexity of actually doing those contrac- 
tions. Naively, one expects that the calculation of the expectation value of a 
quantum state consisting of let us say N spin 1/2's would involve a number of 
operations in the order of the size of the Hilbcrt space and hence exponential 
in N. This strategy would correspond to contracting the tensor network from 
top to bottom. However, we can be smarter and contract the tensor network 
from left to right such that we never have to store more than Di x Di variables 
during the computation (i.e. Di from the upper row and Di from then lower 
row) . Furthermore, the contraction of tensors from left to right can be done in 
a smart way. Let us e.g. assume that all virtual bonds have dimension D and 
physical bonds dimension d] the leftmost contraction corresponds to contracting 
the physical bond, leaving a tensor with two indices i, j of dimension D and D. 
Next we contract one of those with a tensor Ai, and thereafter with Ai. The to- 
tal number of multiplications that had to be done in that procedure is given by 
d 2 D 3 , and this computational cost of contracting the whole network is therefore 
Nd 2 D 3 as we have to repeat the previous step N times. If D is not too big (i.e. 
D < 1000), this can still be done efficiently on present day computers Q. Note 
that it is also obvious how to calculate the expectation value of any observable 
that is defined as a tensor product over local observables [f] ; this is useful to 
calculate e.g. the decay of correlations. 

In the case of periodic boundary conditions, the network to contract is very 
similar (figure [5}:) but one has to keep track of more variables (namely the ones 
going to the left and to the right), and the total computational cost amounts to 
Nd 2 D 5 ; the extra factor of D 2 is responsible for the fact that simulations with 
periodic boundary conditions are considerably slower than with open boundary 
conditions |121[ I137[ I136j . but one has to note that the scaling is still only 
polynomial and that a D in the order of 100 is still achievable in principle, 
which is more than enough to get a good accuracy in actual calculations. 

The following part of this subsection is pretty technical and can be skipped 
for a reader who is not really interested in technical details that are relevant for 
actually implementing numerical renormalization group methods. Let us first 
show how to efficiently calculate expectation values of generic Hamiltonians with 
only terms that act on nearest neighbour spins: 



Instead of calculating the expectation value of each term in the sum separately, 
one can be a bit smarter and store the information of the result of the contraction 
from left to right and arrange everything such that the total computational cost 

Calculations up to D = 10000 have been done in cases where very high precision was 
needed, but such cases are only tractable when explicitely making use of quantum numbers 



8 Note also that any observable can be written as a sum of tensor products; this can easily 
be seen by using e.g. the singular value decomposition 
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is basically still NcPD 3 , which in practice means a huge gain in computing time. 
See Appendix [C] for an actual implementation of this in matlab. 

Another technical tool that is of great practical value is the fact that a matrix 
product state does not have a unique representation in terms of tensors A 1 but 
that so-called gauge transformation on the virtual level leave the physical state 
invariant. This follows from the fact that a multiplication of two matrices is left 
invariant when a nontrivial resolution of the identity is inserted between them: 

AB = AXX~ l B = A'B' 
A' = AX 
B' = X~ X B 

It happens that those gauge degrees of freedom can be exploited to make the 
forthcoming numerical optimization methods better conditioned. In the partic- 
ular case of a MPS with open boundary conditions, we would like them to have 
some orthonormality properties. Consider e.g. again the numerical renormal- 
ization group of Wilson that was sketched in Section [2] There, one obtained 
collections of MPS {IV'cL}} & t each step of the recursion such that they were 
all orthonormal to each other. It is easy to see that a gauge transformation 
can always be found such that this is fulfilled in the case of MPS with open 
boundary conditions. Consider therefore a MPS with tensors A 1 , A 2 , ...A N , and 
let us look at the two collections of Dk MPS defined on the left and right parts 
of the qubits respectively by opening the £/th contracted bonds 

\^n) = A ai A a2 ...A ak .e n \ai)\a 2 ):-\a k ) (7) 

a 1 a 2 ...a k 

\Xn) = A Qk+1 A ak+2 ...A aN \a k+1 )\a k +2)--\a N ) (8) 

cck+ie^cck+2---aN 

where e„ denote the different unit vectors in a D k dimensional vector space. 
Let us furthermore consider the matrices 

A nnl = (VM> (9) 

Bun' = (XM) (10) 

and take their respective square roots A = XX* and B = YY^ . It is clear 
that if we would now do the gauge transformation A\ — ► A\X~ X and repeat 
the procedure described above, we would see that the set {\ipn)} would form an 
orthonormal set as A 1 would be equal to the identity. Similarly, we could make 
the MPS at the right side orthonormal. 

It is interesting to note that the Schmidt coefficients obtained by considering 
the bipartite cut over the fc'th bond are precisely given by the singular values 
of the matrix XY T ; indeed, the matrices A -1 and Y~ T were used to make the 
left and right hand side orthonormal, and the original MPS can of course be 
written as 

71=1 
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Writing out the singular value decomposition of XY T = UYA/ T explicitely, a 
natural gauge transformation to represent the MPS with open boundary condi- 
tions would then be obtained by implementing the gauge-transformation 

A k -> A\X- X U 
A k+1 — > V T Y~ T A k+1 

for all k = 1...7V and by writing out the matrix product state in such a form 
that the Schmidt coefficients are appearing explicitely: 

fofr>= £ AiX 1 AlX 2 ...A» N \a 1 )...\a N ) 

Here the Sj are the diagonal matrices containing the Schmidt coefficients. This 
parametrization of MPS with open boundary conditions coincides exactly with 
the definition of the family of states introduced by G. Vidal in the context of 
the simulation of real time evolution of quantum many-body systems [127] , and 
we have just shown that every MPS with open boundary conditions can written 
like that. 

It would be very appealing to have a similar parametrization in the case of 
MPS with periodic boundary conditions. However, it is clear that no notion of 
singular values or orthonormality can exist in that case as the left side is always 
connected to the right side. However, in a quantum spin system defined on 
a ring with a correlation length much shorter than the length of the ring, the 
boundary effects are expected not to be too pronounced, and an approximate 
orthonormalization can still be carried out [121] . This will be discussed in a later 
section. Related to this issue, let us suppose that we have a MPS with periodic 
boundary conditions for which we know that there exists a MPS description 
with all tensors A 1 equal to each other (the state is hence obviously translational 
invariant). But let us assume that we have a MPS description of that state with 
all tensors A 1 different from each other (i.e. the symmetry has been spoiled 
by site-dependent gauge transformations), e.g. as the result of a variational 
optimization. Is there a way to recover the symmetric description? For this, we 
have to find the gauge transformation A l a — X^A^XZ-^, i = 1 . . .N, such that 

A 1 - A 2 - ■ ■ ■ - A N = A 

Without loss of generality, we can assume that X\ is equal to the identity 
(indeed, we still have the freedom of a global gauge transformation). But then 
we can find X^ and Xjq as the equations A^X^ 1 = X^A^ , a = l..d has at least 
as many equations as unknowns (in practice, this should be done using a least 
squares algorithm). This can then be iterated until all gauge transformation 
Xi have been determined. Alternatively, the new matrices A a can be found by 
considering the gauge-invariant normal form of the MPS: in terms of A l a , the 
norm is expressed as Tr(Ei ■ ■ ■ En), with Ei = X)q=i ® ^L- This expression 
is equal to the norm expressed in terms of A a , namely E , with E = A a ® 
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A a . Therefore, the A a 's are obtained by taking the iV.th root of E\ ■ ■ ■ Em and 
performing a Schmidt-decomposition of the result. 

To summarize, we have shown that MPS are very appealing from the com- 
putational point of view as all correlations functions and the expectation value 
of local Hamiltonians can be calculated efficiently. Moreover, we have shown 
how gauge transformations can be used to bring a MPS with open boundary 
conditions into a normal form. 

A relevant observation is that we can also efficiently calculate the overlap 
between two different MPS. This makes the MPS-approach a very nice tool to 
detect quantum phase transitions, because small changes in the Hamiltonian 
lead to big changes in the ground state around those transition points, and 
those changes can be detected by calculating the overlap between the different 
MPS-approximations of the respective ground states [22] . 

A final remark is that it is also simple to calculate expectation values of 
Hamiltonians with long-range terms. This allows to easily generalize the meth- 
ods introduced in later sections to situations of long-range Hamiltonians. 

3.1.4 Generalization of MPS 

The valence bond state representation can readily be generalized to trees or 
higher dimensions. The generalization to higher dimensions will be discussed 
in Section [6] and there it will become clear that the calculation of expectation 
values is much more involved than in the 1-D case and can only be done approx- 
imately. In contrast to this, it is easy to see that the generalization of MPS to 
tree-networks without loops leads to a family of states from whom all expecta- 
tion values can be calculated efficiently without having to make approximations 



Let us illustrate this for the particular case of a matrix product state on a 
Cay ley tree with coordination number 3 (see figure [6]). To calculate expectation 
values of such a state (^>|0|^>), we start contracting the indices from the bound- 
ary inwards. Because there are no loops, it is easy to see that the number of 
variables to keep track off does not explode exponentially but remains bounded 
just like in the case of matrix product states. Let us for simplicity assume that 
all projectors P are all equal to each other. At any point of the contraction 
when going from the outside to the inside, we have two incoming vertices and 
one outgoing one; as the earlier contractions of the incoming vertices did never 
entangle with each other, their joint state is in a product Pk® Pk- Analogously 
as in the MPS case, the role of the projector P is to apply a completely positive 
map on this product state: if we define the set of Kraus operators A z a .g = P^g^, 



and yields the input for the next level of contractions. As the map from pk 
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Figure 6: A MPS in the form of a Cayley tree with coordination number 3: just 
as in the case of spin chains, the big open circles represent projections of three 
virtual spins to one physical one. 

Pk+i is clearly a nonlinear completely positive map, its structure is much richer 
than in the case of matrix product states and associated completely positive 
linear maps [30] : it happens that such nonlinear maps have multiple fixed points 
with different basins of attractions, and hence the boundary terms can affect 
the expectation values in the bulk of the system (note that this seems to be a 
consequence of the fact that we work with a Cayley tree, where the number of 
vertices increases exponentially as a function of the distance to the center of the 
tree). 

Anyway, the important point to make is that expectation values of any prod- 
uct observable can be readily calculated efficiently by multiplying matrices with 
each other: in the case of a tree with coordination number c, the computational 
cost is D c+1 with D the dimension of the bonds (note that this yields D 3 for 
c = 2 (i.e. MPS) and D 4 for the Cayley tree with c = 3). Note also that it is 
straightforward to obtain a canonical normal form for the tree in the same way 
as we obtained one for the MPS by making use of the singular value decomposi- 
tion [lOlj . Most of the simulation techniques described in the following sections 
can therefore also immediately be generalized to this setting of trees. 

It can be seen that it is also possible to contract the tensor network exactly 
and efficiently if there are loops, but not too many of them. A clear demonstra- 
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tion of this is of course the one-dimensional case of MPS with periodic boundary 
conditions: numbering the d-dimcnsional spins from 1 to 27V in the case of a 
ring with an even number of sites, one can define N new spins by merging spins 
i and 2N — i + 1 into composite pairs. The corresponding MPS with open 
boundary conditions has spin dimension d 2 and bond dimension D 2 , which can 
still be contracted efficiently Similarly, we will be able to contract a generic 
tree with loops efficiently if it can be mapped onto a different tree without 
loops by merging collections of spins together; this procedure will still lead to a 
polynomial computational complexity iff the local dimension of the new spins is 
bounded by a polynomial in the size of the system. This condition is equivalent 
to the requirement that the maximal amount of spins that has been merged into 
a new superspin scales at most logarithmic in the size of the systerro (which 
guarantees a complexity poly(D, TV)). 

One can also formulate a variant of the tree network case, in which the 
nodes of the tree do not carry any physical spins but only the end points of the 
branches do; in other words, we start from a collection of virtual singlets, order 
them in the form of a tree, and put a projector on all the vertices that maps 
the D c dimensional Hubert space to a 1-dimcnsional one. Contraction of such a 
tensor network can again be done in an efficient way by starting the contraction 
at the leaves and working oneself inwards. An interesting observation is the fact 
that every MPS with a given D can be expressed as such a Cayley tree with 
coordination number c = 3 and bond dimension D 2 (this is true both for the 
case of open and periodic boundary conditions 0) . This can easily be seen from 
the construction depicted in figure [7J Such a description was used to develop a 
formalism in which one can do successive renormalization group transformations 
on quantum states as opposed on the usual level of Hamiltonians (X 17] . 

In some further extension, Vidal observed that more general types of tensor 
networks can be contracted efficiently: if all tensors appearing in the tree (with 
only physical spins at the bottom) are isometries, then one can additionally put 
a collection of so-called disentangling unitarics between the different branches. 
The reason why expectation values of local observables (or more generally ob- 
servables that act nontrivially only on a constant number of spins) can still be 
calculated exactly is that in the expression (ip\0\ip) most of those unitaries dis- 
appear cancel each other and hence play no role: one only has to keep track of the 
unitarics within the lightcone of the local observable, and it can be shown that 
the cost of this is still polynomial (albeit with a very large power, especially if 

9 It turns out that this notion of transforming an arbitrary graph into a tree by merging 
vertices together is a well studied problem in graph theory, and that the tree width of a 
graph follows from the optimal way of doing this such as to minimize the maximal number 
of vertices that has to be merged 1951 over all possible mappings from a graph to a tree: the 
tree width is then the maximal number of vertices that is merged in this optimal solution. 
Although calculating the trecwidth of a general graph seems to be NP-hard, efficient (i.e. P) 
approximation algorithms exist that approximate the tree width to within a constant. The 
connection between tree width and MPS was pointed out in |102| . 

10 Note that the opposite situation is different: if we want to represent a generic state on 
the Cayley tree of bond dimension D as a MPS, the MPS will have a site-dependent bond 
dimension bounded above by D l ° 92 ( N ^ = N log2 ^ D ^ (which is still polynomial). 



23 




Figure 7: Every MPS can be represented in the form of a tree network in which 
the vertices are projectors on a 1-dimensional subspace (depicted as square 
boxes as opposed to ellipses for the case of projectors to physical spins); in the 
last step, it is understood that the original matrices A z a p are absorbed into the 
square blocks on the lowest level which are tensors with 6 D-dimcnsional indices. 

one considers the 2-dimcnsional variants of this construction) . This approach is 
called the multi-scale entanglement renormalization ansatz (MERA) |130[ 1128] , 
and has a particularly nice interpretation as successive rescaling transformations 
in the sense of the renormalization group. Due to the presence of the tree net- 
work, the Schmidt number in the 1-D setting can grow logarithmically with the 
system size, and the method is therefore very promising for describing critical 
systems. As opposed to the class of MPS however, there does not seem to be 
a straightforward well conditioned way of doing a variational optimization over 
the class of MERA to find the MERA with minimal energy. Several promis- 
ing approaches have been described in |130[ I128|, 1125] ; the simplest approach is 
to parameterize the unitaries and isometries as exponentials of anti-Hermitean 
operators, and implement the steepest descent optimization algorithm on those 
parameters 0. Further extensions are possible, and there is currently an effort 
in trying to identify the broadest class of quantum states for which all local prop- 
erties can be calculated efficiently (i.e. polynomial complexity). As an example, 

11 Such optimization techniques over the manifold of unitaries have been studied in great 
detail in the context of control theory [?] and are called flow equations; they were also used 
in the context of entanglement theory [4] and first implemented on the level of MERA in 1251 . 
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one can make use of the concept of matrix product operators to show that an 
efficient contraction is still possible when one acts with a quantum circuit of 
arbitrary depth but consisting only of commuting gates (e.g. exp(iaa l z <g> o~ z )) 
between any two spins (even ones that are very far apart from each other) of an 
arbitrary MPS; such states can violate the area law in an extreme sense because 
they might lead to volume laws, and such states might be relevant in simulations 
of systems very far from equilibrium. 

3.2 Reformulating numerical renormalization group meth- 
ods as variational methods in the class of MPS 

Let us now look back at the numerical renormalization group, and reformulate 
it as a variational method within the class of matrix product states. To start 
with, let us aim to find the best possible description of the ground state of e.g. 
the Heisenberg model within the space of all MPSs of the form ([6]), using the 
matrix elements of the matrices {A™} as variational parameters to minimize the 
energy (note that we do not impose the condition that the A n are projectors 
anymore). Using a Lagrange multiplier to ensure normalization, this leads to 
the following optimization problem: 

|^«)G{MPSd} L 

This cost function is multiquadratic in the dN matrices {Ak} with a multi- 
quadratic constraint; indeed, every matrix appears only once in the bra and 
once in the ket, and so this problem is basically equivalent to a cost function of 
the form 

x L ,x z ,... — 

Ki,/e:2,...£i,£2,... 

where the x k are vectors with d x X -D^+i elements and the big tensor Q 
contains the information of the Hamiltonian. The Lagrange constraint is of a 
similar form, and it that case Qk t ,k 2 --- = ^k ± i 2 Sk 2 i 2 ■■■■ There is a standard way of 
solving such an optimization problenT^l which is called alternating least squares 
(ALS). This ALS is an iterative method that works as follows: after making an 
initial guess of the vectors x k , we keep x 2 , ...,x N fixed and optimize over x 1 . 
That subproblcm is of the form 

mmx 1<! HeffX 1 — Xx 1 ^ N e j fx 1 

X 1 

(note that H e ff and N e ff depend on all the other x k ) and is exactly solvable as 
it is a quadratic problem with a quadratic constraints. The solution is that we 
have to choose x equal to the smallest eigenvalue of the generalized eigenvalue 
problem H e ffX = XN e ffX. A crucial point in all this is that there is an efficient 

12 Multiquadratic optimization problems can in principle be NP-hard to solve [7], and so 
there is no guarantee that the alternating least squares method will converge to the global 
optimum. However, in practice this does not seem to occur. See also 1291 . 
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way of calculating H e ff and N e ff given an MPS and the Hamiltonian: these 
are obtained by contracting tensor networks of the type discussed in 13.1.31 and 
an important point in the actual implementation of this is to store the relevant 
matrices for later use. In the next step, we will fix a; 1 , a: 3 , a; 4 , ... and repeat the 
same procedure, until we reach TV and repeat the procedure again by going from 
N to 1 and so on until convergence. Note that at each step of this procedure, 
the energy is going down and we are hence guaranteed to converge to some value 
that is hopefully close to the minimal energy of the Hamiltonian |121j . We will 
call this procedure the variational matrix product state method (VMPS). 

And what about excitations? Within the framework discussed until now, 
this can easily be done variationally [55]. Suppose we found the MPS \ipo) with 
lowest energy, the variational problem is then to find the MPS \ip\) with the 
lowest energy that is orthogonal to the | Vo) - This amounts to taking another 
Lagrange constraint in the optimization problem: 

min [(^il^lVi) - A(^ 1 |^i) - M(^il^o)] 
|V>i)e{MPS D } 

and this can still be solved iteratively in a very similar way. Basically, at every 
step we have a subproblem of the form 

mmx* H e f fx — \x* N e ffX — fiy'x 

whose solution is given by the solution to the generalized eigenvalue problem 
PH e ffPx = XPN e f fPx with P the projector on the subspace orthogonal to 
the vector y. Doing this iteratively, this will again converge to a state \tpi) that 
is the best possible approximation of the first excited state using a MPS. Of 
course, a similar procedure can now be used to construct more excited states. 

Let us now look back at the numerical renormalization group method of Wil- 
son. In essence, the goal is to find the low-energy spectrum of the Hamiltonian. 
NRG can be understood as a different way of doing the optimization discussed 
above, but the NRG method is suboptimal and is only applicable in the par- 
ticular case where there is a clear separation of energies (such as for a Kondo 
impurity or SIAM) . Basically, the steps done in the case of NRG are equivalent 
to the steps one would do using the variational method discussed above during 
the first sweep from left to right. However, the NRG method is stopped after 
that first sweep, which effectively means that one never takes into account the 
influence of the low energy modes on the larger energy ones This is clearly 
suboptimal, and a better result can be obtained by sweeping back and forth 
a few more times until complete convergence is obtained |122j . The computa- 
tional complexity of this is no greater than the one of the original NRG method. 
The method obtained like that is then basically equivalent to the density ma- 
trix renormalization group (DMRG) introduced by S. White [1371 1136] . and it 
is indeed well known that DMRG has a much wider range of applicability than 

13 This feedback may be small in practice, but it is not strictly zero, and its importance 
increases as the logarithmic discretization is refined by taking A — > 1. 
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site index n site index n 

Figure 8: (taken from |122j ) Energy level flow of the S1AM as a function of 
the site index n calculated with (a) VMPS using Z?mps = 32 calculated after 
sweeping and (b) NRG using D N rg = 32 2 = 1024. The blue lines in (a) 
denote the spectrum of <S> 1| + If <S> 7i| ff and demonstrate the decoupling 
of the |- and | spin chains for large n. The inset compares the ground state 
energy as function of D for VMPS (circles) and NRG (squares). The subtracted 
extrapolated ground state energy is E* ~ —3.0714. 

NRG. Note however that NRG and DMRG were always considered to be rather 
different, and it is only by reformulating everything in terms of MPS that the 
many similarities become apparent. 

Let us now come back to the variational method explained above for doing 
NRG. The effective Hamiltonian at chain length n, the central object in NRG, 
is given by H^p = (^Sl^ n |^s)i an d can hence also easily be recovered. Let us 
now illustrate the above by applying them to the SIAM given in equation ([3]). 
The following example is taken out of the paper |122j . Since the Hamiltonian 
couples | and | band electrons only via the impurity, it is possible (see also 
|93j ) to "unfold" the semi- infinite Wilson chain into an infinite one, with f 
band states to the left of the impurity and J, states to the right, and hopping 
amplitudes decreasing in both directions as r-M/ 2 . Since the left and right 
end regions of the chain, which describe the model's low-energy properties, 
are far apart and hence interact only weakly with each other, the effective 
Hamiltonian for these low energies will be of the form 7i | ff <X) 1 ^ + 1 f <X> 7i | ff . This 
is illustrated by the black and blue lines in Fig. [5Ja) . Since 7i| ff and H.^ can 
be calculated separately, instead of simultaneously as in typical NRG programs, 
the dimensions of the effective Hilbert spaces needed in the VMPS approach and 
NRG approaches to capture the low energy properties with the same precision 
are related by Dmps — V f nrg 5 implying significant computational gain with 
VMPS (i.e. a square root speed-up). 

Figure [5] compares the energy level flows calculated using NRG and VMPS. 
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They agree remarkably well, even though we used Dmps = 32 = \/£>nrg- The 
accuracy can still be improved significantly by using larger D and smaller A 
(and exploiting symmetries, not done here for VMPS). 

3.3 Density Matrix Renormalization Group 

As discussed in the previous section, DMRG can effectively be understood as 
the more advanced version of NRG in the sense that one also sweeps back and 
forth. Historically, the discovery of the DMRG algorithm had an enormous 
impact on the field of strongly correlated quantum spin systems as it lead to 
the first algorithm for finding the low energy spectrum of generic spin chains 
with local or quasi- local interactions. One of the first checks of the accuracy 
of the method was the calculation of the ground state energy of the critical 
spin 1/2 Heisenberg chain, and a dazzling precision was obtained when D was 
taken in the order of a few hundreds. The DMRG method also allows one to 
calculate gaps and gave very strong evidence for the existence of the so-called 
Haldane-gap in the case of the spin 1 Heisenberg spin chain. 

Still, the variational matrix product state approach discussed in the previous 
section and the traditional DMRG differ in crucial details, and when looking 
back at the way DMRG has been understood and derived, it is not completely 
obvious to see the parallels with the VMPS approach. We will not explain here 
how DMRG works, but we refer to the many good review papers on the sub- 
ject [57JE0]. O ne important difference is e.g. that the VMPS approach works as 
well in the case of open as of periodic boundary conditions, whereas DMRG uses 
an ansatz that is specifically tailored for systems with open boundary conditions: 
the basic reason for this is that DMRG implicitly uses the orthonormalized nor- 
mal form of MPS discussed in section 13.1.31 but such a normal form does not 
exist in the case of periodic boundary conditions |121j . In practice, simulations 
of spin chains with periodic boundary conditions have been done using DMRG, 
but this is precisely done in the way that was discussed in section 3.1.4 by folding 
the MPS with PBC to one with OBC by squaring the dimension of the virtual 
spins D — ► D 2 ; the cost to pay is an algorithm whose computational cost scales 
as D 6 and, more importantly, for which it is not possible to formulate an ansatz 
of excited states with a definite momentum (see section 3.4). 

Another difference has to do with the fact that the standard DMRG al- 
gorithm is set up in such a way that one performs a variational minimization 
over two sites together instead of over one site described in the VMPS. Such as 
step is then followed by performing a singular value decomposition and throw- 
ing away the least relevant modes; that step is strictly speaking not variational 
anymore, but there is evidence that the conditioning and the numerical conver- 
gence is better when doing it like this. Another reason why this is done is that 
the so-called truncation error gives some idea about the quality of the energy 
estimates obtained in a particular DMRG run. But from the point of view of 
VMPS, the one-site optimization procedure should perform better for the same 
computational cost, as the computational cost in the two-site setup goes up 
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with a factor of d 2 and one could better use that time to increase D. This was 
already observed in the DMRG community Q31 [1351 W\, as DMRG can 
also readily be reformulated with one-site optimizations, but it took until 2005 
and a nice paper S. White [134] before the numerical conditioning problems 
were solved (those techniques of course also apply to VMPS) and the one-site 
DMRG became the standard DMRG method. 

But instead of pointing out the differences, let us discuss a few technical 
improvements of VMPS that parallel the nice features of DMRG. Following 
Fcynman's credo that you should never underestimate the pleasure of readers to 
read what they already know, left repeat again how the whole VMPS procedure 
works. The goal is to minimize 



which is a multiquadratic optimization problem and can be solved using the 
alternating least squares method (ALS). This leads to the problem of solving 
the generalized eigenvalue problem H e ffX = XN e ffX which is a standard opti- 
mization problem itself that is efficiently solvable if the condition number (i.e. 
smallest singular value) of N e ff is not too small. It happens now that in the case 
of open boundary conditions, one can always make appropriate gauge transfor- 
mations (see section l5.1.3[) such as to assure that N e ft is equal to the identity 0. 
Those gauge transformations are always implicitly done in the case of DMRG, 
and should also be done in the VMPS approach as this obviously leads to the 
best conditioning (the computational cost for doing so is small; to see this in 
practice and also to see how easy this VMPS approach is to implement, we refer 
to appendix [Cl for the explicit matlab code for doing so.). 

Let us now see whether there is any good justification for using this ALS- 
method in the context of minimizing the ground state energy of local spin Hamil- 
tonians, both used in the VMPS and in the DMRG approach. Suppose for 
example that we are optimizing over some site in the VMPS method and the 
corresponding D = 128 while the physical spin dimension is d = 2; then we 
are effectively doing an optimization over m = \og d (dD 2 ) + 1 = 15 sites, as 
the number of degrees of freedom we are varying over is equal to the number 
of degrees of freedom in 15 qubits. Ground states of spin models arc typically 
such that there is a finite correlations length (for critical systems, correlations 
are decaying pretty fast too), such that "boundary effects" of more than 7 sites 
away will not play a big role. Colloquially, the bulk convergence to a local op- 
timum cannot occur if the gap in this system of 15 qubits is larger than the 
effect from that boundary, and this gives a handwaving argument why DMRG 
almost always converges to the global minimum. In other words, the success 
of the MPS approach is related to its inherent capability of patching together 

14 Note that in the case of periodic boundary conditions, no such gauge transformation exists 
that makes N e f j equal to the identity, and hence more heuristic methods have to be used. 
We will discuss this in more detail in the section dealing with periodic boundary conditions; 
that makes the VMPS algorithms with periodic boundary conditions a bit more tricky to 
implement, but the basic idea and resulting numerical accuracy is very similar. 




(11) 
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solutions of local (e.g 15 sites) optimization problems, together with the fact 
that MPS are of course rich enough to approximate ground states of arbitrary 
ground states that obey an area law (see appendix B for a proof). 

Before proceeding, let us look at a completely different way of doing varia- 
tional optimizations using MPS when the Hamiltonian under interest is trans- 
lational invariant, and let us consider the thermodynamic limit (i.e. infinitely 
many sites). We know that the ground state will be translational invariant, so 
we might well use an ansatz that reflects this by choosing all A 1 equal to each 
other. This approach has been originally proposed by Rommer and Ostlund [96j 
and later studied by several other authors (see [33] and references therein) . The 
optimization problem is now reduced to minimizing the expectation value of one 
term in the Hamiltonian (i.e. energy per site), and this can be calculated as 
follows: first consider the "transfer matrix" E = A a ®A a and its eigenvalue 
decomposition E = \\ri) (h\ with the A; in decreasing order. The energy 
can now be expressed as 

#=T2 H ai a 2 : a [a> 2 (lo\(A ai ®A a[ )(A a2 (Z>A a , 2 )\ro) 

OL\OLiot^a 2 

where we assumed that the Hamiltonian is only acting on nearest neighbors. 
This cost function is clearly a very nonlinear function of the variables A 1 , and 
standard techniques such as conjugate gradient methods can be used to minimize 
that expression. However, it happens that this optimization procedure may get 
stuck in local minima, and the situation only gets worse when increasing D. In 
comparison with the DMRG or VMPS approach, this method does not seem to 
work very well for several problems. At first sight, this seems to be strange as in 
those latter approaches one has much more variables (i.e. one tensor per site). 
However, this can be understood from the point of view of optimization theory, 
where it is standard practice to introduce more variables than needed such as 
to assure that the problem becomes better behaved. This is precisely what is 
happening here. However, as we will see later, the idea of using translational 
invariant MPS can turned into a successful algorithm by combining it with the 
concept of imaginary time evolution of MPS. It is a bit of a mystery why that 
approach works better than e.g. conjugate gradient methods, but it is probably 
related to the inherent robustness of algorithms evolving in imaginary time. 

As a last remark, we note that it can be very useful from a computational 
point of view to exploit symmetries in the system. This can be done both in 
the DMRG and in the MPS-approach, and we refer to [73J Q3M EM for morc 
details. 

3.3.1 Excitations and spectral functions 

An issue where the usefulness of VMPS really reveals itself is in the context of 
the study of excitations. In the standard DMRG approach, the idea is basically 
to store all the information about the ground state and excited states into 
one big MPS: during the iteration step at e.g. site k, one keeps all tensors 
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A 1 , A 2 , ...A 1 , A k+1 , A N fixed and within the dZ? 2 -dimensional subspace one 
identifies a number of lowest lying states. One then moves to the next site, 
and continues until convergence (again, for a more complete understanding of 
how DMRG works, we refer to the many review papers on the subject [501 157] V 
That procedure is clearly suboptimal and not variational anymore, as there is no 
reason why the same tensors should be used for the ground and excited states. In 
the worst case scenario, the tensors A 1 are block-diagonal, each block containing 
the information of a different excited state, and hence the computational cost 
scales badly with the number of excited states encoded like that (i.e. the cost 
of doing a simulation with k excited states is k 3 times the cost of doing it for 
the ground state; also, the memory requirements scale as k 2 ). 

Looking back at the VMPS approach described above however, this is a 
more natural way to deal with excitations at a lower cost, both with respect to 
memory and computational cost. To repeat again what was explained earlier, 
one can build up the spectrum in a sequential way: first look for the ground 
state, after this start over again and find the MPS with minimal energy that is 
furthermore orthogonal to the ground state, and so further |92| . This procedure 
does not have any of the drawbacks mentioned above in the case of DMRG, and 
is fast and reliable, although it requires to run the whole variational procedure 
for each required excited state. 

If the Hamiltonian under consideration has some symmetries, there might 
of course alternative ways of finding excited states. Consider e.g. the spin 1 
Heisenberg chain. It is well known that the ground state lives in the sector with 
total spin 0, and the really interesting quantity in this case is to find the gap 
between this ground state and the state with minimal energy out of the spin 1 
sector. As the total spin commutes with the Hamiltonian, we can add a small 
magnetic field to the Hamiltonian that plays the role of a chemical potential, and 
within some parameter range, which can easily found by doing some numerics, 
we are guaranteed that the ground state of the complete Hamiltonian will have 
spin 1. It will however still be necessary to implement the procedure mentioned 
above if more excitation energies are to be found. 

Let us continue now and find out whether the formulation of DMRG in terms 
of MPS also allows one to obtain lower bounds to the energies. More precisely, 
suppose the variational method converges to a MPS \ip) with associated energy 
E = (ip\Tl\ip}. Let us next define the quantity 



It is a simple exercise to prove that there exists an exact eigenvalue E ex of TL 
such that E > E ex > E — e. We will show that, in the case of MPS, one can 
calculate the quantity e at essentially the same computational cost as E. This 
implies that the VMPS approach outlined allows to get both upper and lower 
bounds to eigenvalues of a given Hamiltonian H; to our knowledge, this is a 
truly unique feature for a variational method. In typical applications like the 
calculation of the ground state energy of Heisenberg antiferromagnets, e ~ 1CP 8 . 
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Let us now sketch how e can be calculated efficiently |122| . First of all, 
we note that expectation values of tensor products of local observables can be 
calculated by multiplying vectors by matrices: 



Let us now try to evaluate an expression like (tp\ (Tt — E) \ip). For simplic- 
ity, let us assume that Ti represents a spin chains with only nearest neighbor 
couplings. Naively, one expects that one will have to evaluate in the order of 
N 2 expectation values of local observables. There is however a much more ef- 
ficient method: going recursively from the first site to the last one, one keeps 
three D 2 dimensional vectors vo,v\, V2 containing terms with respectively 0, 1, 2 
interaction terms (note that H 2 contains at most two interaction terms). At 
each recursion step, one can easily update vq as a function of vq and the local 
terms of the MPS, and equivalently v\ as a function of v\ and vq plus the local 
Hamiltonian; t>2 can be updated as a function of the local MPS terms, the local 
Hamiltonian and Vo,Vi,V2- Therefore the computational complexity of calcu- 
lating e scales as ND 3 , just as for the case of evaluating the energy {ip\H\^p). 
Another and perhaps more direct way to see this is do make use of the formal- 
ism of matrix product operators (see section 5): there is will be shown that 
every Hamiltonian with only nearest neighbour interaction has a very simple 
parametrization as a matrix product operator, and hence also the square of it. 

As a side product, this insight allows one to devise efficient algorithms for 
calculating highly excited eigenstates and eigenenergies of Hamiltonians. Sup- 
pose for example that one would like to know the closest eigenstate and -energy 
to a certain prespecificd energy E sp . This could be interesting when one wants 
to calculate gaps in the spectrum. The variational problem that has to be solved 
in that case is the following: 



This is indeed again a multiquadratic problem that can be solved using the same 
variational techniques as outlined before, yielding an algorithm with the same 
computational complexity as the original one. 

Similar techniques also allow us to calculate Green's functions in a variational 
way [122] . Green's function have been calculated using NRG methods [2"T1 12"0] . 
but because of the strong interplay between the different energy levels in that 
context, the MPS-approach should be able to give more precise results. The 
DMRG techniques that have been developed for calculating Green's functions 
of spin systems [UJ [571 [S3] are related to our approach but differ in several 
aspects. 

The typical Green's functions of interest are of the form 
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Figure 9: (taken from [122] ) Impurity spectral function for the SIAM, A\ (uj) = 
Ai(lu)ttT/ sin 2 (nd7r/2), normalized such that -4|(0) should be 1, calculated with 
VMPS (solid) and NRG (dashed line). 



where represents the ground state and is a creation operator. Here we will 
assume that \ifj) has been calculated using the MPS-approach, and we showed 
that this can be done to a very good and controlled precision. The basic idea 
is now to calculate the Green's function variationally within the set of MPS by 
finding the unnormalized MPS \x), commonly called a correction vector |122j . 
that minimizes the weighted norm 



\x) 



1 



H 



irj 



w=(n-ui) 2 +r) 2 



Here we used the notation |||£)Hvy = (^I^IOj and the weight W > was intro- 
duced to make the problem tractable (it should not really affect the precision 
as all positive definite norms arc essentially equivalent). Writing \x) in its real 
and imaginary part |x) = \x r ) + i\Xi) an d assuming H, real, this norm can 
be written as 



N = { X r\{H-uf +T 1 2 \Xr)-2( X r\('H-Uj)P\^ 

+( X i\ (H - ujf + rflxt) - 2(Xi\fHi>) + 



(12) 



Minimizing M clearly involves two independent optimizations over |Xr),|Xi)j 
which we will both parameterize as MPS. Both of the optimizations involve 
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minimizing the sum of a multiquadratic and a multilinear term; as terms like 
(x| {H — uj) 2 \x) can be calculated efficiently, we can again use the same tricks 
and keep all but one projectors {P^} fixed and optimize over the remaining one. 
As a multilinear term instead of a quadratic constraint is present, each iteration 
step can be done efficiently by solving a sparse linear set of equations |122j . 
Iterating this procedure until convergence, one can evaluate both \xr), \Xi)i an d 
exactly determine the precision of the result by calculating the norm (fl2|) ; if this 
norm is too large, one can always increase the control parameter D associated 
to the MPS \xr), \Xi)- Finally, one can easily evaluate G(w + irj) in function of 
\x)- The precision of this evaluation can again be bounded. To illustrate this 
with an example, we again consider the SIAM model (see figure [5]). 

3.4 DMRG and periodic boundary conditions 

The variational matrix product state approach discussed in the previous sections 
is both applicable to the case of open (OBC) and periodic boundary conditions 
(PBC). The main difference between both methods is the computational cost 
of the contraction of the associated tensor network (ND 3 versus ND 5 ) and the 
fact that a useful orthonormalization can be used in the case of OBC. On the 
other hand, there are several reasons why simulations with PBC are preferred: 
1. the boundary effects are much smaller than in the case with OBC. This has 
been identified as a serious problem for standard DMRG since its conception, as 
even for very long chains the boundary effects can be seen in the bulk. 2. Due to 
the fact that any nontrivial system with OBC cannot be translational invariant, 
it is not possible to study excitations with a definite momentum. Contrasting 
this to the case of PBC, it is possible there to construct a variational class of 
MPS with a definite momentum, allowing e.g. to obtain a convenient picture of 
the energy-momentum relation. The overall price to be paid however is that the 
computational cost for working with MPS with PBC scales as D 5 as opposed to 
D 3 for OBC. Note that the traditional DMRG algorithm has also been used to 
treat systems with periodic boundary conditions. It is clear, however, from the 
structure of MPS that the only way to represent a MPS with PBC by one with 
OBC is to use the virtual bonds to "teleport" a maximally entangled state from 
the first to the last site. This leads to an effective D that is equal to the original 
one squared (this is the only possibly way to accomodate the extra degrees of 
freedom to take into account this teleportation), and hence the computational 
cost for the same accuracy would scale as D 6 . But even in doing so, it would 
be very hard to construct states with a definite momentum. 

When implementing the VMPS approach to minimize the energy of a given 
Hamiltonian with periodic boundary conditions |121| . special precautions have 
to be taken such that the problem remains well conditioned. In practice, this 
means that we have to make sure that the matrix N e ff in the generalized 
eigenvalue problem H e ffX = XN e ffX is positive definite and that its smallest 
eigenvalue is as big as possible. There are several options to achieve this. First 
of all, we can still play with the gauge conditions which would correspond to a 
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Figure 10: (taken from |121j ) Comparison between DMRG (squares) [136] and 
the (new) VMPS method (circles) for PBC, and N = 28. For reference the 
DMRG results [136] for the Heisenberg chain with OBC (triangles) arc also 
shown. Inset: variation of the local bond strength from the average along the 
chain, calculated with the new method and D = 40. 



transformation 

N eff (X <g> Y <g> I d ) N eff (X* <g> Y f <g> I d ) . 

Looking back at the case with OBC, the reason why we could always make 
N e ff — I is that the left and right side of the chain factorize such that N e ff 
was always a tensor product to start with. Here, the PBC enforce the existence 
of correlations between those sides, hence enforcing N e ff to be "correlated". 
However, in practice the amount of correlations will be small (as these are finite 
size effects), and hence N e ff will be close to a tensor product. A sensible 
way of choosing the gauge transformation therefore consists out of two steps: 
1. find the tensor product N\ (g> N% Cg> Id that approximates best N e ff (this 
can easily be achieved by doing a singular value decomposition in the space 
of Hcrmitcan operators); 2. choose the gauge transformations X,Y such that 
XNiX^ = I = YN2Y'. During the first sweeps, when absolute accuracy is not 
relevant yet and N e f / could still be far from a tensor product, one could e.g. add 
some identity to N e fj such as to make it better conditioned (N e ff — > N e ff+el). 
The procedure outlined here is only one of many possible ones, and it might well 
be that another choice of gauge conditions leads to better performance, but the 
one described here seems to give good results in practice |121j . An example is 
provided in figure [TOl where we looked at the Heisenberg spin 1/2 chain with 
N = 28 (we choose this example with a small number of sites as we can still 
compare this with exact results, and for reference also the results of a DMRG 
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Figure 11: (taken from [92]) Lowest states of a bilinear-biquadratic 5 = 1 
chain, N = 40 sites, D = 10. (a) 9 = -tt/2, E a = -2.7976N, (b) 9 = -0.74tt, 
Eq = — 1.4673^. Empty circles: lowest energy states. Filled circles: first branch 
of excitations. We estimate an absolute error AEk ~ 10~ 3 , by comparison with 
calculations with larger D. 

calculation for the same Hamiltonian is provided). As seen in the figure, the 
scaling of the error for MPS with PBC nicely follows the scaling of the error 
obtained in a DMRG-calculation for a Hamiltonian with OBC. 

As already explained in the previous section, it is easy to look for excited 
states using an extra Lagrange constraint. However, the VMPS approach can 
also easily be generalized such as to become a varational method over states with 
a definite momentum. The basic idea is very simple. Consider the state |92j 

N-l ikn 
n— v 

| X ) = Y,Tr{A} tl A* a ...A» N )\ ai )\ a2 )...\ aN ) 

aii... 

where the operator Tfe is the shift operator implementing a translation over k 
sites. The state \ipk) has the momentum k by construction and is obviously a 
superposition of N MPS (note that this implies that the block entropy of the 
state \ipk) can scale as \og(ND 2 ) as opposed to log(D 2 ) for a normal MPS). 
What is clear, however, is that the whole VMPS procedure outlined above can 
now be repeated for this extended class of MPS: the energy is still a multi- 
quadratic function of all tensors A 1 involved, and expectation values of local 
observables can still be calculated efficiently by contracting the corresponding 
tensor network. The price to pay is an extra factor of N (the number of sites) 
in all the calculations, due to the fact that cross terms between the different 
superpositions enter the optimization. To achieve this rather mild slowdown, 
a rather involved scheme of bookkeeping has to be maintained, and this is ex- 
plained in great detail in the paper [92] . As an illustration of the power of this 
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technique, figure [TT] represents the calculation of the energy-momentum relation 
for two values of the bilinear-biquadratic S = 1 spin chain parameterized by the 
Hamiltonian 

H = J2 cos(6)SiSi + i + sm(9) (SiS i+1 
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4 Time evolution using MPS 



One of the big advantages of the formulation of renormalization group methods 
in terms of matrix product states is that it also allows to describe time evolution. 
This opens up the possibility of simulating the non-equilibrium properties of spin 
chains, a topic that is currently very relevant given the recent breakthroughs of 
creating strongly correlated quantum spin systems in e.g. optical lattices [T7] , 
Since the development of DMRG, several methods have been proposed [HI I7T1 
[13 H3H1 HHl [2H HIH I3H E3 [H] . In the spirit of this paper, we will only review 
the variational approach, and refer to the review of Schollwock and White for 
DMRG-related approaches [98] . 

4.1 Variational formulation of time evolution with MPS 

Mathematically, the problem is to evolve an initial MPS in real time by updating 
the tensors in the MPS-description under Hamiltonian evolution. In practice, 
this could be used to investigate how excitations travel through the spin chain or 
to get spectral information about a Hamiltonian of interest. In a similar spirit 
as the previous sections, we would like to do this in a variational way: given a 
Hamiltonian and an initial MPS, evolve that state within the manifold of MPS 
in such a way that the error in approximating the exact evolution is minimized 
at every infinitesimal step [119j . 

We are particularly interested in the case where the Hamiltonian, which can 
be time-dependent, is a sum of local terms of the form 



and where < ij > means that the sum has to be taken over all pairs of nearest 
neighbours (note that the situation with long-range interactions can be treated 
in a very similar way). There are several tools to discretize the correspond- 
ing evolution. This is not completely trivial because generically the different 
terms in the Hamiltonian don't commute. A standard tool is to use the Trotter 
decomposition [1091 1110] 



Suppose e.g. that the Hamiltonian can be split into two parts A and B such that 
all terms within A and within B arc commuting: H = A + B: A = J^i At; B = 
Bi] [Ai, = = [Bi, Bj}_. This ensures that one can efficiently represent 
e tA as a product of terms. The evolution can then be approximated by evolving 
first under the operator e lStA , then under e iStB , again under e lStA and so further. 
The time step can be choosen such as to ensure that the error made due to this 
discretization is smaller than a prespecified error, and there is an extensive 
literature of how to improve on this by using e.g. the higher order Trotter 
decompositions |105[ [82] . In the case of nearest neighbour Hamiltonians, a 
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a) A= H H H H 



B= 



I 1 = cxp(i(T. H> (T ) 

b) A= 



Figure 12: Spin network representation of the different operators obtained in 
the case of a Trotter expansion of the Ising Hamiltonian in transverse field 
with a) nearest neighbour decomposition and b) Ising versus magnetic field 
decomposition. 

convenient choice for A is to take all terms that couple the even sites with the 
odd ones to the right of it and for B the ones to the left of it. In that case, 
e lA and e %B are tensor products of nearest-neigbour 2-body operators; see figure 
112a for a pictorial representation of this in terms of spin networks. However, 
it is important to note that different choices of A and B might work better 
practice. Consider e.g. the Ising Hamiltonian in an transversal field and its 
decomposition into two different terms A and B (see also Fig [12b ) : 

n Ising = ® + H E °x ■ ( 13 ) 



Obviously, a similar kind of decomposition is possible in the case of the Heisen- 
berg model, but there three terms are needed instead. The advantage of evolu- 
tion with such a decomposition is that it does not break translational invariance 
as in the even-odd case. 

Let us next investigate how to treat this time-evolution in a variational way 
within the class of MPS. A sensible cost function to minimize is given by 

\\A\m)-m+i))\\i (14) 

where \ip(k)) is the initial MPS, \i/j(k + 1)) is the one to be found, and A is the 
operator arising out of the Trotter expansion 0. First of all, it is important 
to note that this cost function can indeed be evaluated efficiently when \ip(h)) 



15 Remark that as alternative cost-function we could have chosen ^ ^^^VJi^lt 1 ^?! which 
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and \ip(k + 1)) arc MPS and A is of the form discussed above. In principle, this 
cost function can be made equal to zero by increasing the bond dimension D 
by a factor of at most d 2 (this is indeed the largest possible Schmidt number 
of an operator acting on two sites). However, the whole point of using MPS is 
that MPS with low bond dimension are able to capture the physics needed to 
describe the low-energy sector of the Hilbert space. So if we stay within that 
sector, the hope is that the bond dimension will not have to be multiplied by a 
constant factor at each step (which would lead to an exponential computational 
cost as a function of time), but will hopefully saturate. This is certainly the case 
if we evolve using e.g. imaginary time evolution of a constant local Hamiltonian, 
as we know that in that case the ground state is indeed well represented with 
a MPS with not too large bond dimension. It is however important to keep 
in mind that an exponential explosion is in principle possible for other kinds 
of evolution: the worst-case computational complexity for time evolution using 
MPS is exponential as a function of time. 

Taking this into account, a justified way of dealing with time evolution is to 
prespecify an error e that can be tolerated, and then look for the minimal D for 
which there exists a MPS \ip(k + l)) that yields an error smaller than e. Looking 
back at the cost function (fT"4"]) , it looks pretty familiar how to minimize it: the 
cost function has only quadratic and linear terms, and wc will hence be able to 
minimize it by a method very similar to the alternating least squares method 
discussed in the previous section. More specifically, the cost function has only 
multiquadratic and multilinear terms in the variables of the MPS, and we will 
solve this in a recursive way where at each recursion step an optimization of the 
form 

x^AeffX - 2x t y eff 

has to be solved. The solution to this problem is the simple solution of the 
linear set of equations 

Aefjx = y eff , 

and this hence leads to a very efficient way of minimizing the cost function 
for time evolution in a recursive way: sweeping back and forth, we solve the 
above optimization subproblem at each step, and we are guaranteed that the 
total cost function goes down at every step and hence will converge. After 
convergence, we can check how big the error has been for the particular value 
of D that we choose, and if this error is too big, we can increase D and repeat 
the optimization. Note that the only requirement for the complete Trotter step 
evolution to be successful is that the tensor network when sandwiching the 
operators A and B between two MPS can be contracted efficiently. It is clear 
that this method both works in real and imaginary time evolution, and both 
for time-independent and time-dependent Hamiltonians. As a note, and this 
will be crucial in the applications of this variational time evolution of MPS to 
higher dimensional spin systems, we also want to point out that nowhere we 

has the aim at finding the normalized MPS \ip(k + 1)} that has maximal overlap with the 
evolved version of the original one. It is however an easy exercise to find out that this 
optimization problem leads exactly to the same optimal MPS, up to a normalization factor. 
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used the fact that the operator A was close to the identity; i.e. this method is 
applicable in a more general way than only for small Trotter steps. Note also 
that this method is applicable to systems with both open and periodic boundary 
conditions, and also with long-range interactions. 

4.1.1 time-evolving block-decimation 

There has recently been a lot of attention on studying time evolution using 
MPS, but instead of using the optimal variational way described above, the 
vast majority of the works has been using the so-called time-evolving block- 
decimation (TEBD) procedure introduced by Vidal |127j . The reason for this is 
that this method can readily be implemented with existing DMRG code [1391 
124] . It explicitely uses the Schmidt normal form described in section l3.1.31 and 
is hence more suitable for MPS with open boundary conditions. Also, one has 
to use the Trotter expansion with even-odd sites decomposition. It can also be 
understood as a variational method, but in the more restricted sense in which 
we consider the following set-up: given a MPS and one (possibly non-unitary) 
gate acting on two nearest- neighbors, find the new MPS with given dimension 
D that approximates this optimally. This can be done using the singular value 
decompisition: using the normal form for MPS with OBC, we know that the 
original MPS is of the form 

IV>) = X>£>lx£> 

n=l 

with {iV'n} an d {IV'n} orthogonal sets. The gate which acts on the nearest 
neighbours locally increased the dimension of the bond with a factor of at most 
d 2 , and we would like to reduce the dimension of that bond to D again. One 
can orthonormalizc everything again and obtain the Schmidt decomposition over 
that particular bond, and then the reduction can trivially be done by discarding 
(i.e. projecting out) the smallest Schmidt coefficient^. In the next step, evo- 
lution between the next nearest neighbours is considered, and so further. Note 
that the variational method discussed above could deal with all those gates at 
once; nevertheless, the computational cost of both methods is very similar, and 
in practice both methods seem to achieve a similar accuracy if the time steps 
are small. The first one, however, can be improved by choosing longer time 
steps and higher Trotter orders [37] . 

It is also possible to devise more sophisticated methods that somehow make 
use of the fact that superpositions of MPS can still be dealt with in an efficient 
way. For a nice review that compares all such methods and the ones described 
before, see [37]. For nice examples of the power of real-time evolution, we refer 
to the review article [98] . 

16 This is indeed the virtue of the singular value decomposition (SVD): if we want to approx- 
imate a given matrix with one of lower rank in an optimal way (in the sense as to minimize 
the Frobenius norm of the difference), then the solution is given by taking the SVD and put 
the smallest singular values equal to zero. 
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4.2 Finding ground states by imaginary time evolution 



The tools discussed in the previous section apply both to real and imaginary 
time evolution. This provides a completely different way of finding ground states 
of local quantum spin Hamiltonians: if we start with an arbitrary state IV'o}) 
evolution in imaginary time leads to 



with H = J2k=i ^k\Xk)(Xk\ the eigenvalue decomposition of H. Hence, as long- 
as the ground state is not degenerate and for times longer than the inverse gap, 
the state \ip(t)) will converge exponentially fast to the ground state, and the 
speed of convergence is exactly quantified by the gap. This is indeed something 
that is a recurring theme: the smaller the gap, the slower all variational methods 
seem to converge. This is of course not unexpected because it is difficult to 
discriminate excited states with small energy from the ground state. However, 
it is interesting to note that the closer a system is to criticality (i.e. gap goes 
to zero), the bigger D has to be for a MPS to approximate the ground state 
for a fixed error (see appendix B to see that in critical systems, this scaling is 
polynomial in the number of spins in the worst case scenario). A very interesting 
question would be to relate the density of states above the gap to the decay of 
the Schmidt coefficients that one gets by dividing the ground state into two 
pieces. 

In practice, finding a ground state using imaginary time evolution is pretty 
reliable. Of course, the time steps have to be adjusted such that they become 
smaller and smaller, but one of the great features of imaginary time evolution 
is its inherent robustness: it does not really matter if one makes mistakes in the 
beginning, as there is anyway an exponential convergence to the ground state 
afterwards. This is in contrast to the evolution in real time. 

4.2.1 Infinite spin chains 

Another advantage of the imaginary time evolution approach is that one can 
easily treat systems with periodic boundary conditions or treat the thermody- 
namic limit (number of sites — > oo) by making use of the inherent translational 
invariance of MPS with all tensors A 1 equal to each other. 

To illustrate how this can be done, let us consider the specific case treated in 
[30] and assume that we want to model the ground state of the Ising Hamiltonian 
in a transverse field defined on a spin chain. As discussed in section [4~T1 following 
equation[]21 one can choose the decomposition in the Trotter step in such a way 
that both the operators A and B 



n 



k=l 



A 



fc+i 



A- 



B = 
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are completely translational invariant, and hence we can stay within the mani- 
fold of translational invariant states to describe its evolution. The next step is 
to see how exp(StA) and exp(5tB) look like. The latter is simple, as it is just 
equal to 

exp(dtB) = <E)kCxp(Stha x ). 

The former expression is obviously a product of commuting operators, and a 
simple exercise allows one to see that there is a simple matrix product description 
for it: 

exp(StA) = Tr(C ll C l2 ...)X n ®X l2 ® ... (1 

Co = (n °a) C t - 



a 








a 


1 








1 



Xo = " x i 



^sinh((5i) cosh(^) 

v /sinh((5t) cosh(St) 

1 
-1 

This can justifiably be called a matrix product operator (MPO) [see next sec- 
tion]. The associated bond dimension is 2, and this means that when acting 
on a MPS of dimension D with this MPO, the exact representation of the new 
MPS will at most be 2D. In this particular example of the Ising Hamiltonian, 
imaginary time evolution would now amount to act subsequently with exp(StA) 
and exp(StB) on an initial state. Combined together, it happens that their 
product Q = exp(StB /2) cxp(StA) exp(5tB/2) is again exactly of the form (|15|1 
with unchanged tensor Ci and Xi but where we have to replace Xo = I with 
Xq = cxp(Sta x ). So Q is a very simple MPO with bond dimension 2, and the 
goal is to evolve a translational invariant MPS using this MPO. This can be 
done is a very simple way: given a MPS with tensor Ai of dimension D, the 
action of Q is such that we have to replace 

A^Y. Ak ® Cl ^\ Xl \ k )- 

k,l 

The new MPS corresponding to this Ai has bond dimension 2D, and this has 
to be reduced as otherwise its size would increase exponentially. The optimal 
choice is again simple if we consider a system with open boundary conditions: 
we consider its normal form (section 13.1.31) . and cut it in the appropriate way 
without losing translational invariancel 17 !. That last step is not so easy to justify 



17 Moro specifically, a possible procedure is as follows I3UI : calculate the leading left \vi) 
and right eigenvector \v r ) of J^ fe Cj, ® Cfc (note that this can be done in a sparse way). 
Reshape those eigenvectors in the square matrices X\ and X r , and as those matrices are the 
fixed points of the CP-maps J2 fe Ck-Cj. and J^ fe Cl.C^, Xi and X r are guaranteed to be 
positive. Next take the singular value decomposition of \/ X r \JXi = t/SW, and reduce the 
number of columns of U and V to D and discard the D lowest singular values in S (U and 
V hence become 2D X D matrices and S a D X D matrix) . Define G r = \/ X^ 1 U VT, and 
G; = V^Wv/jfj -1 , and make the transformation Ai — > G[.AiG r such that it corresponds to 
a MPS with D-dimcnsional bonds instead of 2D. 
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rigorously, but seems to work very well in practice: what we do is calculating the 
Schmidt normal form with respect to the 2D dimensional bonds, and then cut all 
bonds together. The tricky thing is that cutting the bond somewhere changes 
the Schmidt coefficients somewhere else, but the point is that these changes 
are only of the order of the Schmidt coefficients that are cut and those are very 
small anyway. Amazingly, this procedure works very well, and even a small bond 
dimension of D = 32 already reproduces the ground state energy-density E of 
the critical Ising model (h = 1) up to a precision better than E — 4/ir < 10~ 7 . 
Much better conditioning can also be obtained |30| by working with hermitean 
matrices {A 1 }. 

Clearly, this procedure can be repeated for any translational invariant Hamil- 
tonian like the Heisenberg model (note that the bond dimension will be bigger 
there), and the big advantage is that it allows to treat infinite systems. The 
finite system with periodic boundary conditions can be dealt with in a similar 
way, although the cutting procedure is more involved there. 

A variant of this procedure can be obtained by using the even-odd Trotter 
decomposition. In that case, exact translational invariance is broken, but it is 
still perfectly translational invariant with period 2: 

\i>) = ^2 ■■■A il B i2 A i3 B ii ...\ii}\i 2 )\i a )\H) 

This type of imaginary time evolution has been studied in detail by G. Vidal 
|129j . and convergence seems to be fast. The updating and cutting works in a 
very similar way as discussed in the example discussed above. Note that from 
the point of view of variational states, it could be advantagous to work with 
such an ABAB... scheme when an antiferromagnetic ordering is expected. 
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5 Matrix Product Operators 



Instead of restricting our attention to pure matrix product states, we can readily 
generalize the approach and deal with matrix product operators (MPO). In its 
most general case, a MPO is defined as [1 191 1144] 

6= ^{A\Al..)a n ® ai2 ®... (16) 

ill2 . . . 

with Oi a complete single particle basis (e.g. the Pauli matrices for a qubit). We 
already encountered an example of such a MPO in the previous section, where 
the evolution following a Trotter step was expressed by a MPO; that MPO 
essentially played the role of a transfer matrix during the evolution. As we 
will show later, any transfer matrix arising in the context of classical partition 
function will have an exact representation in terms of such a MPO. This is 
the reason why all the rcnormalization methods described above can also be 
used in the context of 2-D classical spin systems. As a side remark, it is also 
true that any translational invariant spin Hamiltonian with nearest neighbour 
interactions has an exact MPO representational. Matrix product operators 
have been shown to be very useful to obtain spectral information about a given 
Hamiltonian. Examples are the paramctrization of Gibbs states (section l5.1|) . 
the simulation of random quantum spin systems (section 15. 2p . the calculation of 
classical partition functions (section 15. 3p , and the determination of the density 
of states (section [5T]) . 

If a MPO is a positive operator and has trace one, then it becomes a matrix 
product density operator (MPDO). A simple systematic way of constructing 
MPDO [119j is by making use of the fact that every mixed state can be seen 
as a part of a bigger pure system, namely its purification. If we model this 
purification as a MPS, then the MPDO is obtained by tracing over the purifying 
degrees of freedom. This picture is especially useful if the purification is such 
that to every original spin, a locally accompanying purifying spin is present, and 

18 For this, let us consider the spin 1/2 case. We first need the fact that there always exists 
a basis such that the Hamiltonian is of the form 

a,i j 
where O can be any one-qubit operator. It is now a small exercise to prove that 

tl«2 — 

Xo = I Xi = a x X2 = cr y X3 = <j z X4 = O 

vi = |0> v r = |4> 

A = |0><0| + |4>(4| Ai = |0)(l| + |l)(4| A 2 = |0>(2| + |2><4| 

A 3 = |0><3| + |3>(4| A 4 = |0><4| 

Actually, one can easily prove that D = 5 is optimal because this is the operator Schmidt 
number of the Hamiltonian when splitting it into two pieces. 
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that the pure state is such that these pairs of spins correspond to one bigger 
site: 



E ^( A in A t 

i>\: 



In>|ji)®l*2)|j2>®... 



(17) 
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This representation of matrix product density operators in terms of purifications 
will turn out to be very useful to describe Gibbs states in thermal equilibrium. 
Note that not all MPO can arise from locally tracing out an ancilla: the ten- 
sors E l arising in the above description are very special in the sense that they 
correspond to completely positive maps and hence live in a convex positive cone. 

5.1 Finite temperature systems 

Matrix Product Operators are very convenient to approximate Gibbs states of 
local Hamiltonians. More specifically, we would like to approximate the operator 



-/3H 



13 n, Un, 



where / is the identity operator. In the spirit of last sections, this amounts to 
evolving the maximally mixed state in imaginary time for a period (3/2 |119[ 
1144] . However, instead of implementing this evolution on the MPO directly, a 
square root speed up in computational complexity can be gained^! by consid- 
ering the evolution on its purification (see previous section). More specifically, 
the initial state / can be seen as a collection of halve of maximally entangled 
pairs. In the notation of equation (fT?]) . the inititial state is a MPS with D = 1 
and where Aij = 8ij . We can now use the variational algorithm discussed above 
to evolve this initial state over a time (3/2, and once we have that state it can 
be used to calculate e.g. the free energy and any correlation function. 

Note that without time evolution, it would be very hard to do a variational 
calculation over all possible matrix product operators to approximate the Gibbs 
state, as that state is variationally characterized by the condition that it mini- 
mizes the free energy Tr (pH) — TS(p); the problem here is the calculation of the 

19 Note that such a square root speed-up is guaranteed when the temperature is very low as 
there as the number of classical correlations is zero for pure states. 
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von-Neumann entropy S(p), which we do not know how to calculate for matrix 
product density operators in general. But as a by-product of the imaginary time 
approach, we can readily calculate the entropy of the Gibbs state p 

Tr (He- 1311 ) 

S( P ) =lo g Tr(e-e»)+P ^ H) K 

as all the quantities on the right hand side of this equation can efficiently be 
calculated and approximated using the MPS/MPO representation obtained via 
evolution in imaginary time. 

Let us now see whether it is possible to calculate non-equilibrium properties 
of such MPO by evolving them. The most general type of evolution allowed 
by quantum mechanics is described by the Lindblad equation, which takes into 
account noise and can therefore increase the entropy. In general, this means that 
the system is interacting with an extra auxiliary one, and this poses somehow 
problems for treating the evolution via the purification described above: the 
dimension of the ancilla has to grow continuously. The straightforward remedy 
is to work with the full matrix product operator picture, and evolution can again 
be treated in a straightforward variational way [1191 H44] . 



5.2 Random quantum spin systems 

Another nontrivial application of matrix product operators and more specifically 
of the corresponding purifications is the fact that those allow to simulate random 
quantum spin systems in a highly efficient way. By making effective use of the 
entanglement present between system and an appropriately chosen ancilla, it 
happens that one can simulate exponentially many realizations of the random 
system in one run; this is quantum parallelism in the strongest possible sense. 
To explain this, we closely follow the exposition given in the paper [86j . 

Let us consider a quantum system with Hilbert space 7i that evolves ac- 
cordingly to a Hamiltonian H(r\, . . . , r n ) where r\, . . . ,r n are random variables 
that take values within a finite discrete set, rg € = {X[ , . . . A^ l£ }, with a 
probability distribution given by p{r\, . . . , r„). In order to simulate exactly the 
dynamics of such a system one would need to perform 11™= 1 m £ simulations, one 
per each possible realization of the set of random variables r = (ri, . . .r n ). For 
each realization the system evolves to a different state \tp r (t)) = e~ lff W i ' /ft |V>o}j 
where \ipo) is the initial state. Given this set of evolved states and a physical 
observable O, one is typically interested in the average of the expectation values 
of that observable in the different evolved states, that is, in quantities of the 
form: 

((d(t))):=J2p(*)(Mt)\o\Mt))- (18) 

r 

On the following we describe an algorithm that allows to simulate in parallel all 
possible time evolutions of the random system described above. We consider an 
auxiliary system with Hilbert space TL a and a Hamiltonian acting on 7i ® Ti. a 
of the form H = H{R\ 1 . . . , i? n ), where R\, . . . , R n are operators that act in 
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7i a , commute with each other and have spectra Ti, . . . , r„. Note that we have 
replaced the set of random variables r by a set of quantum operators R with the 
same spectra. The algorithm works as follows. 1) Initialization. Let us prepare 
the auxiliary system in an initial superposition state of the form: 

li)=Ev^)l r ). ( 19 ) 

r 

where the states |r) are simultaneous eigenstates of the set of operators R, with 
i?£|r) = ri\v). Each state |r) is therefore in one to one correspondence with one 
realization of the the set of random variables r, its weight in the superposition 
state (|19[) being equal to the probability with which the corresponding realiza- 
tion occurs for the random system. 2) Evolution. We evolve the initial state of 
the composite system |-!/>o) <8> IV'a) under the Hamiltonian H. The evolved state 
is 

!*(*)> =EVpWllM*)>®|r>. (20) 

r 

This superposition state contains the complete set of evolved states we are 
interested in. 3) Read-out. In order to obtain the quantities (|18j) we just need 
to measure the observable O ® 1, 

(*(i)|d®i|*(i)> = (<d(t)>>. (2i) 

The algorithm above allows us, in particular, to obtain the averaged properties 
of a random system over the collection of all possible ground states. Let us 
assume that the interaction between the system and the ancilla is introduced 
adiabatically, so that the Hamiltonian is now H(t) = _ff(/3(i)R), where (3(t) is 
a slowly varying function of time with /3(0) = 0, (i{T) = 1, T being the time 
duration of the evolution. If the system is prepared in the ground state of the 
Hamiltonian H(0), the algorithm above will simulate in parallel all possible 
adiabatic paths, so that the composite superposition state (|20[) will contain all 
possible ground states of the random system [86j . 

Additionally, the scheme above can be easily extended for the computation 
of other moments of the distribution of physical observables (higher than (fl~8|) ). 
which are sometimes important in the understanding of QRS [86] , For example, 
quantities like ((O 2 ) - (O) 2 ), can be computed by using an additional copy of 
the system. 

The algorithm described above reduces the simulation of a quantum random 
system to the simulation of an equivalent non-random interacting problem. This 
exact mapping allows us to integrate the simulation of randomness in quantum 
systems within the framework of numerical methods that are able to efficiently 
simulate the corresponding interacting problem. As an illustrative example we 
consider the case of a ID spin s = 1/2 system with random local magnetic field. 
The Hamiltonian of the system is: 

N 

H{h,...,b N ) = H + Bj2kS!, (22) 
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Figure 13: Numerical simulation ol the time evolution of a random field XY spin 
chain with N = 40. We show the correlation function ((c\.Ck)) as a function 
of time and momentum k. Here ct oc sin(k£)cg, k = j^ri , • ■ ■ , jfxj an( i 
Cf = Yl^i Sg,(Sf + iS v t ) are the fermionic operators given by the Jordan- 
Wigner transformation. The evolution Hamiltonian is H(bi, . . . , 6jy) = Ho + 
B^2^ =1 bgSf with H being the XY Hamiltonian with B = 0, B/J = —4 and 
p(h) = 1/2 N . The initial state |V>o) is the ground state of Hq. As the system 
evolves in time the initially sharp Fermi sea disappears. 
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where Hq is a short range interaction Hamiltonian, b = (61, . . . , bjy) is a set 
of classical random variables that take values {1/2, —1/2} with probability dis- 
tribution p(h). Following the algorithm above the 2 N simulations required for 
the exact simulation of the dynamics (or the ground-state properties) of this 
random problem can be simulated in parallel as follows. We consider an aux- 
iliary ID spin cr=l/2 system. We prepare this ancilla in the initial state 
\ip a ) = X)b a b|k), where the states |b) have all z components of the N spins 
well defined, <rf |b) = bi\h), and «b = \Jp(h). The entangled properties of the 
state of the ancilla reflect the classical correlations among the random variables. 
For example, for a uniform distribution of the random field, p(h) = l/2 N , the 
state of the ancilla is just a product state, \ij) a ) oc (| f) + | 1))® N ■ We evolve the 
system and the ancilla under the interaction Hamiltonian 

H = Hffi, ...,a z n )=H + (3j2 Sf 51. (23) 

t 

Here, f3 = B if we want to simulate dynamics under Hamiltonian (|22|) . and 
(3 is a slowly varying function of time with (3(0) = and (3(T) = B for the 
simulation of the ground state properties. We have then reduced the simulation 
of the random problem to that of the time evolution of two coupled spin 1/2 
chains with Hamiltonian (|23|) . This problem is equivalent to a ID lattice problem 
of N sites with physical dimension d = 2x2, which can be easily incorporated to 
the framework of the variational numerical methods introduced in the previous 
sections. 

5.3 Classical partition functions and thermal quantum states 

In this section, we will show how the concept of MPO can be used to calcu- 
late the free energy of a classical 2-dimensional spin system. This also leads to 
an alternative way of treating thermal states of 1-D quantum spin systems, as 
those can be mapped to classical 2-D spin systems by the standard classical- 
quantum mapping r°l Historically, Nishino was the first one to pioneer the use 
of DMRG-techniques in the context of calculating partition functions of classi- 
cal spin systems [78]. By making use of the Suzuki-Trotter decomposition, his 
method has then subsequently been used to calculate the free energy of trans- 
lational invariant 1-D quantum systems [111 11031 1132] , but the main restriction 
of those methods is that it cannot be applied in situations in which the number 
of particles is finite and/or the system is not homogeneous; furthermore, one 
has to explicitely use a system with periodic boundary condition in the quan- 
tum case, a task that is not well suited for standard DMRG. The variational 
MPS-approach gives an easy solution to those problems. 

2t) In the context of MPS and especially PEPS, there also exist a different mapping between 
classical and quantum spin models in the same dimension 11241 . There, the thermal classical 
fluctuations map onto quantum ground state fluctuations, and this leads to a lot of insights 
in the nature of quantum spin systems. 
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Our method relies in reexpressing the partition and correlation functions as 
a contraction of a collection of 4-index tensors, which are disposed according to 
a 2-D configuration [74] , We will perform this task for both 2-D classical and 

1- D quantum systems. 

Let us consider first the partition function of an inhomogencous classical 

2- D n-lcvcl spin system on a L\ x Li lattice. For simplicity we will concentrate 
on a square lattice and nearest-neighbor interactions, although our method can 
be easily extended to other short-range situations. We have 

Z = J2 exp[-pH(x n ,...,x L ^)], 

X 11 ,...,X L l L 2 

where 

H (x 11 , ...) = [Hf (x ij ,x i+1 ' j ) + H% (x l \x^ +1 ) 

ij 

is the Hamiltonian, x lJ = 1, . . . , n and (3 is the inverse temperature. The singular 
value decomposition allows us to write 

n 

exp [-/3HV(x,y)] =Y.fU x )aUy)> 

a=l 

with q G {4, — >}. Defining the tensors 

n 

2 = 1 

the partition function can now be calculated by contracting all 4-index ten- 
sors X lJ arranged on a square lattice in such a way that, e.g., the indices 
l,r,u,d of X^ are contracted with the indices r,l,d,u of the respective ten- 
sors P -1 ^, X t+1 ' 3 . In order to determine the expectation value 
of a general operator of the form 0({x l: >}) = ZW^ O 1 ^ (x 1 ^ ) , one just has to 
replace each tensor X^ by 

n 

*iL «>") = E^w^w^ 1 "'^)/-^)^ 1 - 

x=l 



The quantum case if very similar. We consider the partition function of an 
inhomogencous 1-D quantum system composed of L n-level systems, 

Z = trexp (-/?#). 

It is always possible to write the Hamiltonian H as a sum H = J^ fe Hk with 
each part consisting of a sum of commuting terms. Let us, for simplicity, assume 
that H = Hi + H2 and that only local and 2-body nearest neighbor interactions 
occur, i.e. Hk = Y^i an d <^' I+1 , 0^ +1 = 0, with i,j = 1,...,L. 
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The more general case can be treated in a similar way. Let us now consider a 
decomposition 



exp 



M° k 



E s k 



fi+i 

ka 



(24) 



The singular value decomposition guarantees the existence of such an expression 
with k < n 2 . As we will see later, a smart choice of H = J2k ^ k can typically 
decrease k drastically. Making use of the Suzuki-Trotter formula Ell 



Tr 



I cxp 




1 

17 



it can be readily seen that the partition function can again be calculated by 
contracting a collection of 4-indcx tensors X 1 ^ defined as 



^(W){rr')ud — 1 ll D lr 1 2l' D 2r 



■id] 



where the indices (1,1') and (r,r') arc combined to yield a single index that 
may assume values ranging from 1 to k 2 . Note that now the tensors X % i and 
X % i coincide, and that the indices u of the first and d of the last row have 
to be contracted with each other as well, which corresponds to a classical spin 
system with periodic boundary conditions in the vertical direction. A general 
expectation value of an operator of the form O = ZO 1 <%> ■ ■ ■ <g> O n can also be 
rccxprcssed as a contraction of tensors with the same structure: it is merely 
required to replace each tensor X 1 ^ in the first row by 



X (U')(rr')ud (^) 



'.id] 



Let us next move on to explain how the tensor contraction can be done. We 
will use the techniques that were originally developed in the context of PEPS 
in order to contract the tensors X lJ introduced above in a controlled way. The 
main idea is to express the objects resulting from the contraction of tensors along 
the first and last column in the 2-D configuration as matrix product states and 
those obtained along the columns 2, 3, . . . , L — 1 as matrix product operators. 
More precisely, we define 



'X 1 



X 1 



E tr(Xll...X^){ ri ...r M \ 
n...rM=l 

rn 

£ tr{^...xZ L )\h...l M ) 



7J 



h,ri 



21 Note that in practice, it will be desirable to use the higher order versions of the Trotter 
decomposition. 
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where m = n for 2-D classical systems and m = k 2 for 1-D quantum systems. 
These MPS and MPOs are associated to a chain of M m-dimensional systems 
and their virtual dimension amounts to D = n. Note that for 2-D classical 
systems the first and last matrices under the trace in the MPS and MPO reduce 
to vectors. The partition function (and similarly other correlation functions) 
reads Z = (X 1 | X 2 • ■ ■ | X L ). Evaluating this expression itcratively by 
calculating step by step (X- 7 | := (X- 7 ^ 1 | X- 7 for j = 2, . . . , L — 1 fails because 
the virtual dimension of the MPS (X 5 | increases exponentially with j. A way 
to circumvent this problem is to replace in each iterative step the MPS (X- 7 | by 

a MPS ( X" 7 | with a reduced virtual dimension D that approximates the state 

(X- 7 | best in the sense that the norm 5K := 1 1 ( X J " | — (X 3 ||| is minimized. Due 
to the fact that this cost function is multiquadratic in the variables of the MPS, 
this minimization can be carried out very efficiently; the exponential increase 
of the virtual dimension can hence be prevented and the iterative evaluation 
of Z becomes tractable, such that an approximation to the partition function 

can be obtained from Z ~ ( X | X L ) . The accuracy of this approximation 
depends only on the choice of the reduced dimension D and the approximation 
becomes exact for D > D L . As the norm SK can be calculated at each step, 
D can be increased dynamically if the obtained accuracy is not large enough. 
In the worst case scenario, such as in the NP-complete Ising spin glasses [5], 
D will probably have to grow exponentially in L for a fixed precision of the 
partition function. But in less pathological cases it seems that D only has 
to grow polynomially in L; indeed, the success of the methods developed by 
Nishino |78j in the translational invariant case indicate that even a constant D 
will produce very reliable results. 

We will illustrate this with an example of bosons in optical lattices. A system 
of trapped bosonic particles in a 1-D optical lattice of L sites is described by 
the Bose-Hubbard Hamiltonian [53] 

L-l jj L L 

H = - J ^(ajdj+i + /i.e.) + hi(hi - 1) + Vifii, 

i—l i—1 i—1 

where a\ and are the creation and annihilation operators on site i and 
hi = ajflj is the number operator. This Hamiltonian describes the interplay 
between the kinetic energy due to the next-neighbor hopping with amplitude J 
and the repulsive on-site interaction U of the particles. The last term in the 
Hamiltonian models the harmonic confinement of magnitude V% = Vo(i — io) 2 . 
The variation of the ratio U/J drives a phase-transition between the Mott- 
insulating and the superfluid phase, characterized by localized and delocalized 
particles respectively |36j . Experimentally, the variation of U / J can be realized 
by tuning the depth of the optical lattice [S31 QH] • On the other hand, one typ- 
ically measures directly the momentum distribution by letting the atomic gas 
expand and then measuring the density distribution. Thus, we will be mainly 
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interested here in the (quasi)-momentum distribution 



1 L 



r.s — 1 



Our goal is now to study with our numerical method the finite-temperature 
properties of this system for different ratios U/J. We thereby assume that 
the system is in a thermal state corresponding to a grand canonical ensemble 
with chemical potential [i, such that the partition function is obtained as Z = 
^ Te -P{H-tiN) _ jjcre, ]\[ — ^\ =1 hi represents the total number of particles. 
For the numerical study, we assume a maximal particle-number q per lattice 
site, such that we can project the Hamiltonian H on the subspacc spanned by 
Fock-states with particle-numbers per site ranging from to q. The projected 
Hamiltonian H then describes a chain of L spins, with each spin acting on a 
Hilbert-space of dimension n = q + 1. A Trotter decomposition that turned out 
to be advantageous for this case is 



i 



M 2 



(25) 



with H = H R +Hs+H T , H R = - + Zti R {l) R {l+1 \ H s = -+ ^ti S^S^ l+1 \ 
Ht = Eti T {1 \ JJ« = at + ~a h = -i(aj - d t ), T« = in,(n, - 1) + V.h, 
and V = e~^ Hrt e~™ Hs e~™( HT ~^\ a\, di and fii thereby denote the pro- 
jections of the creation, the annihilation and the number operators aj, and m 
on the g-particle subspace. The decomposition (|M|) of all two-particle operators 
in expression (|25[) then straightforwardly leads to a set of 4-index tensors X^ ud , 
with indices I and r ranging from 1 to (q + l) 3 and indices u and d ranging 
from 1 to q + 1 . Note that the typical second order Trotter decomposition with 
H = H cvcn + -ff dd would make the indices I and r range from 1 to (q + l) 6 . 

Let us start out by considering the limit U/J — » oo in which double occu- 
pation of single lattice sites is prevented and the particles in the lattice form a 
Tonks-Girardeau gas [57]. In this limit, the Bose-Hubbard Hamiltonian maps 
to the Hamiltonian of the exactly solvable (inhomogeneous) XX-model, which 
allows to benchmark our algorithm. The comparison of our numerical results 
to the exact results can be gathered from fig. Q3] Here, the density and the 
(quasi)-momcntum distribution are considered for the special case f3J = 1, 
L = 40, N = 21 and Vq/J~ 0.034. The numerical results shown have been 
obtained for Trotter- number M — 10 and two different reduced virtual dimen- 
sions D = 2 and D = 8. The norm SK was of order 10 -4 for D = 2 and 10~ 6 
for D = 8 [3- From the insets, it can be gathered that the error of the numer- 
ical calculations is already very small for D = 2 (of order 10 -3 ) and decreases 
significantly for D = 8. This error can be decreased further by increasing the 
Trotter- number M. 



22 We note that we have stopped our iterative algorithm at the point the variation of SK 
was less than 10 -8 . 
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Figure 14: (taken from |74j ) Density and (quasi)-momentum distribution in 
the Tonks-Girardeau gas limit, plotted for (3 J = 1, L = 40, N = 21 and 
Vq/J = 0.034. The dots (crosses) represent the numerical results for D = 2 
(D = 8) and the solid line illustrates the exact results. From the insets, the 
error of the numerical results can be gathered. 

As the ratio U/J becomes finite, the system becomes physically more inter- 
esting, but lacks an exact mathematical solution. In order to judge the reliability 
of our numerical solutions in this case, we check the convergence with respect 
to the free parameters of our algorithm (q, D and M). As an illustration, the 
convergence with respect to the parameter q is shown in figure 1151 In this fig- 
ure, the density and the (quasi)-momentum distribution are plotted for q = 2, 3 
and 4. We thereby assume that (3 J = 1, L = 40 and N — 21 and consider 
interaction strengths U/J = 4 and 8. The harmonic potential Vb is chosen in a 
way to describe fib-atoms in a harmonic trap of frequency Hz (along the lines 
of |87j). We note that we have taken into account that changes of the ratio U /J 
are obtained from changes in both the on-site interaction U and the hopping 
amplitude J due to variations of the depth of the optical lattice. The numerical 
calculations have been performed with M = 10 and D = q + 1. From the figure 
it can be gathered that convergence with respect to q is achieved for q > 3. 
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Figure 15: (taken from [74] ) Density and (quasi)-momentum distributions for 
interaction strengths U/J = 4 and 8. Here, (3 J — 1, L — 40, N — 21 and 
M = 10. Numerical results were obtained for q = 2 (plus-signs), q = 3 (crosses) 
and q = 4 (solid line). For comparison, the distributions for U/J = (dotted 
lines) and U/J — > oo (dash-dotted lines) are also included. 



5.4 Density of States 

Information about the density of states can be obtained by studying the quantity 

f(t) = Tr (e~ im ) 

as a function of time. Indeed, assume that we know the function f(t) exactly 
from t = — oo — > +oo. Its Fourier transform 

with {Ek} the spectrum of the Hamiltonian H . Hence the Fourier transform of 
fit) gives all the information about the spectrum. 

Of course, in practice we will not be able to determine fit) for all times, and 
we will only be able to approximate it within some time window. The effect 
on its Fourier transform is that it becomes convolved with a sinc-function of 
some width inversely proportional to the time window This means that the 

23 Of course, we can also make use of more advanced signal processing tricks to get a better 
conditioning. We refer to 1831 for more details. 
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resolution of such simulations will depend on the time frames in which we can 
calculate f(t). Furthermore, only a discrete number of points will be available, 
such that we have to do a discrete Fourier transform instead, but this is still 
sufficient to get the spectral information we're interested in. 

As first discussed by T. Osborne, the calculation of f(t) can efficiently be 
implemented with the techniques discussed before [53]. Basically, it is the real- 
time equivalent of calculating e~@ H , and as shown in the previous sections, 
there are several options for doing this. First of all, we can evolve a collection of 
maximally entangled states with 7 <£> e~ lHt using small Trotter steps, and then 
calculate its overlap with the original maximally entangled states (note again 
that several options exist for how to choose the Trotter decomposition). Indeed, 
it happens that 

(I\A®I\I) = TrA 

with 1 7) = J2k 1^)1^) a maximally entangled state. 

Alternatively, we can consider the Suzuki- Trotter decomposition for an in- 
finitesimal step, write down the complete tensor network corresponding to f(t), 
and contract it using the techniques discussed in the previous section f5. 31 Note 
that in this particular case, there is again the possibility of contracting the 
corresponding tensor network in two directions, namely in the time or space 
direction. Experience will have to tell us which way is more reliable. 

In any case, it is really exciting to see how many new tools and possibilities 
can be explored by reformulating it in the MPS-language. 
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6 Projected Entangled Pair States and ground 
states of 2-D quantum spin systems 



In this section, we present a natural generalization of the ID MPS to two and 
higher dimensions and build simulation techniques based on those states which 
effectively extend DMRG to higher dimensions. We call those states projected 
entangled-pair states (PEPS) [1181 [75] . since they can be understood in terms 
of pairs of maximally entangled states of some auxiliary systems that are locally 
projected in some low-dimensional subspaces. This class of states includes the 
generalizations of the 2D AKLT-states known as tensor product states [5TJ EU 
[79l [80l \77\ 172"! 1114] which have been used for 2D problems but is much broader 
since every state can be represented as a PEPS (as long as the dimension of 
the entangled pairs is large enough). We also develop an efficient algorithm to 
calculate correlation functions of these PEPS, and which allows us to extend the 
ID algorithms to higher dimensions. This leads to many interesting applications, 
such as scalable variational methods for finding ground or thermal states of spin 
systems in higher dimensions as well as to simulate their time-evolution. For the 
sake of simplicity, we will restrict to a square lattice in 2D. The generalization 
to higher dimensions and other geometries is straightforward. 

We want to emphasize however that the PEPS method is not yet as well 
established as the 1-D MPS or DMRG methods; this is mainly due to its bigger 
complexity, but also to a big extent due to the fact that these PEPS methods 
are relatively unexplored which makes that there is a lot of way for improvement 
and exciting research. 

6.1 Construction and calculus of PEPS 

There have been various attempts at using the ideas developed in the context 
of the numerical renormalization group and DMRG to simulate 2-D quantum 
spin systems. However, in hindsight it is clear why those methods were never 
very successful: they can be reformulated as variational methods within the 
class of 1-dimcnsional matrix product states, and the structure of those MPS 
is certainly not well suited at describing ground states of 2-D quantum spin 
systems. This can immediately be understood when reconsidering the area law 
discussed in the first section (see figure |T6|) : if we look at the number of degrees 
of freedom needed to describe the relevant modes in a block of spins, this has to 
scale as the boundary of the block, and hence this increases exponentially with 
the size of that boundary. This means that it is impossible to use a NRG or 
DMRG approacrF^. where the degrees of freedom is bounded to D. 

However, it is straightforward to generalize the MPS-picture to higher di- 
mensions: the main reason of the success of the MPS approach is that it allows 

24 Clearly, the VMPS/DMRG methods can reveal very valuable information in the case 
of quasi 2-dimensional systems such as ladders with a few rungs; see e.g. 11381 for a nice 
illustration. 
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Figure 16: Representation of a quantum spin system in 2 dimensions using 
the PEPS representation. If we calculate the entropy of a block of spins, then 
this quantity will obey an area law and scale as the length of the boundary 
between the block and the rest. PEPS-states are constructed such as to have 
this property build in. 

to represent very well local properties that are compatible with e.g. the transla- 
tional symmetry in the system. These strong local correlations are obtained by 
sharing maximally entangled states between neighbours, and the longer range 
correlations are basically mediated by the intermediate particles. This is of 
course a very physical picture, as the Hamiltonian does not force any long- 
range correlations to exist a priori, and those only come into existence because 
of frustration effects. This generalization to higher dimensions can therefore be 
obtained by distributing virtual maximally entangled states between all neigh- 
bouring sites [114j . and as such a generalization of the AKLT-picturc is obtained. 

More specifically, each physical system at site i is represented by four aux- 
iliary systems at, fej, Cj, and di of dimension D (except at the borders of the 
lattices). Each of those systems is in a maximally entangled state 

D 

|/>=£|m> 

i=l 
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Figure 17: Structure of the coefficient related to the state | fcn, .., &44 ) in the 
PEPS I^a)- The bonds represent the indices of the tensors [Ai] k that are 
contracted. 

with one of its neighbors, as shown in the figure. The PEPS | "J ) is then obtained 
by applying to each site one operator Qi that maps the four auxiliary systems 
onto one physical system of dimension d. This leads to a state with coefficients 
that are contractions of tensors according to a certain scheme. Each of the 
tensors is related to one operator Qi according to 

[Mwud = {k\Qi\l,r,u,d) 

and thus associated with one lattice site i. All tensors possess one physical 
index k of dimension d and four virtual indices Z, r, u and d of dimension D. The 
scheme according to which these tensors are contracted mimics the underlying 
lattice structure: the four virtual indices of the tensors are related to the left, 
right, upper and lower bond emanating from the corresponding lattice site. The 
coefficients of the PEPS are then formed by joining the tensors in such a way 
that all virtual indices related to same bonds are contracted. This is illustrated 
in fig. [T7]for the special case of a 4 x 4 square lattice. Assuming this contraction 
of tensors is performed by the function J r (-), the resulting PEPS can be written 
as 

d 

ki,...,k M = l 

This construction can be generalized to any lattice shape and dimension and 
one can show that any state can be written as a PEPS if we allow the bond 
dimension to become very large. In this way, we also resolve the problem of 
the entropy of blocks mentioned above, since now this entropy is proportional 
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to the bonds that connect such block with the rest, and therefore to the area of 
the block. Note also that, in analogy to the MPS [33], the PEPS are guaranteed 
to be ground states of local Hamiltonians. 

There has recently been a lot of progress in justifying this PEPS picture; 
M. Hastings has shown [35] that indeed every ground state of a local quantum 
spin Hamiltonian has an efficient representation in terms of a PEPS, i.e. one 
whose bond dimension D scales subexponentially with the number of spins 
under interest. Also, he has shown that all thermal states have an efficient 
representation in terms of matrix product operators. This is great news, as it 
basically shows that we have identified the relevant manifold describing the low- 
energy physics of quantum spin systems. This can lead to many applications 
in theoretical condensed matter physics, as the questions about the possibility 
of some exotic phase of matter can now be answered by looking at the set of 
PEPS hence skipping the bottleneck of simulation of ground states. 

The family of PEPS also seems to be very relevant in the field of quantum 
information theory. For example, all quantum error-correcting codes such as 
Kitacv's toric code [59j exhibiting topological quantum order have a very simple 
and exact description in terms of PEPS [124j . Furthermore, the PEPS-picturc 
has been used to show the equivalence between different models of quantum 
computation [1 14j ; more specifically, the so-called cluster states [M] have a 
simple interpretation in terms of PEPS, and this picture demystifies the inner 
workings of the one-way quantum computer. 

6.2 calculus of PEPS 

We now show how to determine expectation values of operators in the state | VP ). 
We consider a general operator O = J7, 0« an d define the D 2 x D 2 x D 2 x D 2 - 
tensors 

[E?]%\^= E (^\n[A*]t ud [A 3 ]l, u , dl . 

k,k' = l 

In this definition, the symbols (W), (rr r ), (uu') and (dd 1 ) indicate composite 
indices. We may interpret the 4 indices of this tensor as being related to the 4 
bonds emanating from site j in the lattice. Then, (4'|0|4') is formed by 
joining all tensors Ej 3 in such a way that all indices related to same bonds are 
contracted - as in the case of the coefficients of PEPS. These contractions have 
a rectangular structure, as depicted in fig. [TgJ In terms of the function J-(-), 
the expectation value reads 

(* |0|tt) =F(E°\...,E° N ). 

The contraction of all tensors £7 • 3 according to this scheme requires a number 
of steps that scales exponentially with N - and makes calculations intractable 
as the system grows larger. Because of this, an approximate method has to be 
used to calculate expectation values. 
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Figure 18: Structure of the contractions in (^a \ ^a)- In this scheme, the 
first and last rows can be interpreted as MPS | Ui ) and ( U4 | and the rows in 
between as MPO U2 and U3. The contraction of all tensors is then equal to 
(KilCWalE/i). 

The approximate method suggested in [118j is based on matrix product 
states (MPS) and matrix product operators (MPO). The main idea is to inter- 
pret the first and last row in the contraction-scheme as MPS and the rows in 
between as MPO. The horizontal indices thereby form the virtual indices and 
the vertical indices are the physical indices. Thus, the MPS and MPO have 
both virtual dimension and physical dimension equal to D 2 . Explicitly written, 
the MPS read 

I EM = E tr([<"] : 

d 1 ,...,d L = l 

(Ul\ = £ tr([E°?] 
fii,...,«i=i 

and the MPO at row r is 

Ur= £ tr([2#*] 

ui,...,u L = l 
d 1 ,...,d L = l 

In terms of these MPS and MPO, the expectation value is a product of MPO 
and MPS: 

(* |0|* ) = (U L \U L - 1 ---U 2 \U 1 ) 



..■[E^] Ul -)\d u ...,d L ) 
■..[E° L ^r 1 ){u l ,...,u L \ 



■■■[E^f LdL )\u 1 ,...,u L )(d 1 ,...J L \. 
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The evaluation of this expression is, of course, intractable. With each multi- 
plication of a MPO with a MPS, the virtual dimension increases by a factor 
of D 2 . Thus, after L multiplications, the virtual dimension is D 2L - which is 
exponential in the number of rows. The expression, however, reminds of the 
time-evolution of a MPS. There, each multiplication with a MPO corresponds 
to one evolution step. The problem of the exponential increase of the virtual 
dimension is circumvented by restricting the evolution to the subspace of MPS 
with a certain virtual dimension D. This means that after each evolution step 
the resulting MPS is approximated by the "nearest" MPS with virtual dimen- 
sion D. This approximation can be done efficiently, as shown in [119] . In this 
way, also ( \0\ ) can be calculated efficiently: first, the MPS | U2 ) is formed 
by multiplying the MPS \Ui) with MPO U 2 . The MPS \U 2 ) is then approx- 
imated by I U2 ) with virtual dimension D. In this fashion the procedure is 
continued until | Ul-i ) is obtained. The expectation value ( & \0\ ) is then 
simply 

(*|0|*} = (t/i I £/l_i>. 

Interestingly enough, this method to calculate expectation values can be 
adopted to develop very efficient algorithms to determine the ground states of 
2D Hamiltonians and the time evolution of PEPS by extenting DMRG and the 
time evolution schemes to 2D. 



6.3 Variational method with PEPS 

Let us start with an algorithm to determine the ground state of a Hamiltonian 
with short range interactions on a square Lx L lattice. The goal is to determine 
the PEPS I W ) with a given dimension D which minimizes the energy: 

(H) = [ - 1 1 ; (26) 

Following |121j , the idea is to iteratively optimize the tensors A4 one by one while 
fixing all the other ones until convergence is reached. The crucial observation is 
the fact that the exact energy of | 3/ ) (and also its normalization) is a quadratic 
function of the components of the tensor A4 associated with one lattice site i. 
Because of this, the optimal parameters Ai can simply be found by solving a 
generalized eigenvalue problem. 

The challenge that remains is to calculate the matrix-pair for which the 
generalized eigenvalues and eigenvectors shall be obtained. In principle, this is 
done by contracting all indices in the expressions ( \& \ H\ \& ) and ( $ | $ ) except 
those connecting to Ai. By interpreting the tensor Ai as a <iD 4 -dimensional 
vector Ai, these expressions can be written as 

(*|ff|*> = AH, A, (27) 
= A\MAi. (28) 
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Thus, the minimum of the energy is attained by the generalized eigenvector Ai 
of the matrix-pair (Hi,Afi) to the minimal eigenvalue /i: 

H t A t = ptMiAi 

It turns out that the matrix-pair {Hi, Mi) can be efficiently evaluated by 
the method developed for the calculation of expectation values: Mi relies on the 
contraction of all but one tensors Ej (with / denoting the identity) according to 
the same rectangular scheme as before. The one tensor that has to be omitted 
is E\ - the tensor related to site i. Assuming this contraction is performed by 
the function Gi(-), Mi can be written as 

ML I"* = 

If we join the indices (klrud) and (k'l'r'u'd'), we obtain the dD A x cLD 4 -matrix 
that fulfills equation To evaluate Gi(-) efficiently, we proceed in the same 

way as before by interpreting the rows in the contraction-structure as MPS and 
MPO. First, we join all rows that lie above site i by multiplying the topmost 
MPS | Ui) with subjacent MPO and reducing the dimension after each multi- 
plication to D. Then, wc join all rows lying below i by multiplying ( Ul | with 
adjacent MPO and reducing the dimension as well. We end up with two MPS 
of virtual dimension D - which we can contract efficiently with all but one of 
the tensors Ej lying in the row of site i. 

The effective Hamiltonian Hi can be determined in an analogous way, but 
here the procedure has to be repeated for every term in the Hamiltonian (i.e. in 
the order of 2N times in the case of nearest neighbor interactions). Assuming 
a single term in the Hamiltonian has the tensor-product structure H s = JT^ h\, 
the effective Hamiltonian H\ corresponding to this term is obtained as 

The complete effective Hamiltonian Hi that fulfills equation (|T7| is then pro- 
duced as 

Hi = ^W; s . 

s 

Thus, both the matrices Mi and Hi are directly related to the expressions 
Qi{E[, ...fEjf} and Gi(E^ 1 , E^™). These expressions, however, can be eval- 
uated efficiently using the approximate method introduced before for the cal- 
culation of expectation values. Therefore, the optimal Ai can be determined, 
and one can proceed with the following site, iterating the procedure until con- 
vergence. 

6.4 Time evolution with PEPS 

Let us next move to describe how a time-evolution can be simulated on a PEPS. 
We will assume that the Hamiltonian only couples nearest neighbors, although 
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more general settings can be considered. The principle of simulating a time- 
evolution step is as follows: first, a PEPS | ^f^) with physical dimension d = 2 
and virtual dimension D is chosen as a starting state. This state is evolved by 
the time-evolution operator U = e~ lHSt (we assume ft = 1) to yield another 
PEPS I * B ) with a virtual dimension Db increased by a factor r\: 

\*b) = U\*° a ) 

The virtual dimension of this state is then reduced to D by calculating a new 
PEPS I ) with virtual dimension D that has minimal distance to \^b)- This 
new PEPS is the starting state for the next time-evolution step. The crucial 
point in simulating a time-evolution with PEPS is thus the development of an 
efficient algorithm for reducing the virtual dimension of a PEPS. 

Before formulating this algorithm, let us recite how to express the product 
U\ ^ A ) in terms of a PEPS. This is done by means of a Trotter-approximation: 
first, the interaction-terms in H arc classified in horizontal and vertical accord- 
ing to their orientation and in even and odd depending on whether the interac- 
tion is between even-odd or odd-even rows (or columns). The Hamiltonian can 
then be decomposed into a horizontal-even, a horizontal-odd, a vertical-even 
and a vertical-odd part: 

H = H he + H ho + H ve + H vo 

The single-particle operators of the Hamiltonian can simply be incorporated 
in one of the four parts (note that different Trotter decompositions are again 
possible, e.g. grouping all Pauli operators of the same kind in 3 different groups 
as we discussed earlier, and in some cases this leads to a clear computational 
advantage). Using the Trotter-approximation, the time-evolution operator U 
can be written as a product of four evolution-operators: 

jj = p -iHSt ^ e -iH h<! 8t e -iH ho 5t e -iH vc 8t e -iH vo St ^9) 

Since each of the four parts of the Hamiltonian consists of a sum of commut- 
ing terms, each evolution-operator equals a product of two-particle operators 
Wij acting on neighboring sites i and j. These two-particle operators have a 
Schmidt-decomposition consisting of, say, 77 terms: 

v 

p=l 

One such two-particle operator applied to the PEPS | ^ A ) modifies the 
tensors A® and A® associated with sites i and j as follows: assuming the sites i 
and j are horizontal neighbors, A® has to be replaced by 

d 

[-^*] l(rp)ud = X/ \- U i\ k> ] Irud 
fe'=l 
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and A° becomes 

k d k 1 

k' = l 

These new tensors have a joint index related to the bond between sites i and j. 
This joint index is composed of the original index of dimension D and the index p 
of dimension r/ that enumerates the terms in the Schmidt-decomposition. Thus, 
the effect of the two-particle operator Wij is to increase the virtual dimension 
of the bond between sites i and j by a factor of 77. Consequently, e - lH ^t anc j 
e -iH ho st mcrease t ne dimension of every second horizontal bond by a factor of 77; 
e -iH V e8t an( j e -iH va st j-Jq ^jjg same f or every second vertical bond. By applying 
all four evolution-operators consecutively, we have found an approximate form 
of the time-evolution operator U that - when applied to a PEPS | ^>° A ) - yields 
another PEPS | b ) with a virtual dimension multiplied by a constant factor 77. 

The aim of the approximate algorithm is now to optimize the tensors A{ 
related to a PEPS | $ a ) with virtual dimension D, such that the distance 
between | f ^) and | ^s) tends to a minimum. The function to be minimized 
is thus 

K(A U ...,A M ) = |||*a) -|*B>f- 
This function is non-convex with respect to all parameters {Ai, Am}- How- 
ever, due to the special structure of PEPS, it is quadratic in the parameters Ai 
associated with one lattice site i. Because of this, the optimal parameters A; 
can simply be found by solving a system of linear equations. The concept of 
the algorithm is to do this one-site optimization site-by-site until convergence 
is reached. 

The coefficient matrix and the inhomogeneity of the linear equations system 
can be calculated efficiently using the method developed for the calculation of 
expectation values. In principle, they are obtained by contracting all indices in 
the expressions for the scalar-products ( ^/a | ^>a ) and ( ^>a | *b ) except those 
connecting to A^, By interpreting the tensor Ai as a cLD 4 -dimensional vector 
Ai, these scalar-products can be written as 

(VaI^a) = A\M t A t (30) 

(*a|*b) = A\Wi. (31) 

Since 

K = ( * B I * B ) + ( *A I *A ) - 2Re{ *A I *B }, 
the minimum is attained as 

MiAi = IV,. 

The efficient calculation of Mi has already been described in the previous 
section. The scalar product ( ^ a | ^ b ) and the inhomogeneity W, are calculated 
in an efficient way following the same ideas. First, the DDb X DDb X DDb X 
DZ?B-tensors 

d 

\ f I ( uu ')i dd ') _ r a*i k rn i fe 

L i\ (W)(rr>) ~ Z^i I J) Irud i n V r' u' d' 
k=l 
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are defined. The scalar-product (^a\^b) is then obtained by contracting 
all tensors Fj according to the previous scheme - which is performed by the 
function 

(y A \y B }=F(F 1 ,...,F M ) 

The inhomogenity Wi relies on the contraction of all but one of the tensors Fj , 
namely the function Gi('), in the sense that 

D 

I' r' u' d' = 1 

Joining all indices (klrud) in the resulting tensor leads to the vector of length 
dD A that fulfills equation (f5T|) . Thus, both the scalar-product ( ^>a \^b) and 
the inhomogenity Wi are directly related to the expressions T{F\, Fm) and 
Qi[Fi, Fm)- These expressions, however, can be evaluated efficiently using 
the approximate method from before. 

Even though the principle of simulating a time-evolution step has been re- 
cited now, the implementation in this form is numerically expensive. This is 
why we append some notes about how to make the simulation more efficient: 

1. - Partitioning of the evolution: The number of required numerical operations 
decreases significantly as one time-evolution step is partitioned into 4 substcps: 
first the state | ^ A ) is evolved by e - %H vo$t only and the dimension of the in- 
creased bonds is reduced back to D. Next, evolutions according to e ~ lH "^ St ^ 
e -iH ho st an( j e -iH hB 5t f ii ow Even though the partitioning increases the num- 
ber of evolution steps by a factor of 4, the number of multiplications in one 
evolution step decreases by a factor of r] 3 . 

2. - Optimization of the contraction order: Most critical for the efficiency of the 
numerical simulation is the order in which the contractions are performed. We 
have optimized the order in such a way that the scaling of the number of mul- 
tiplications with the virtual dimension D is minimal. For this, we assume that 
the dimension D that tunes the accuracy of the approximate calculation of Mi 
and Wi is proportional to D 2 , i.e. D = k.D 2 . The number of required multipli- 
cations is then of ordciF^I k 2 D w L 2 and the required memory scales as dr]K 2 D 8 . 

3. - Optimization of the starting state: The number of sweeps required to reach 
convergence depends on the choice of the starting state for the optimization. 
The idea for finding a good starting state is to reduce the bonds with increased 
virtual dimension r/D by means of a Schmidt-decomposition. This is done as 
follows: assuming the bond is between the horizontal neighboring sites i and j, 
the contraction of the tensors associated with these sites, Bi and Bj, along the 
bond i—j forms the tensor 

Dr) 

[ Ml j\lud r'u'd' = X] [ Bi \lpud[ B j]pr'u'd'- 
p=l 

25 The scaling D 10 is obtained when at all steps in the algorithm, a sparse matrix algorithm 
is used. In particular, we have to use an iterative sparse method for solving the linear set of 
equations in the approximation step. 
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By joining the indices (klud) and (k'r'u'd'), this tensor can be interpreted as a 
dD 3 x cLD 3 -matrix. The Schmidt-decomposition of this matrix is 



dD i 

P =\ 

with the Schmidt-coefficients c p (c p > 0) and corresponding matrices A^ and 
Aj. We can relate these matrices to a new pair of tensors A® and A® associated 
with sites i and j: 



1 Ipud V P J lud 



U0] K = 

L J J prua V r L j J rud 

The virtual dimension of these new tensors related to the bond between sites i 
and j is equal to the number of terms in the Schmidt-decomposition. Since 
these terms are weighted with the Schmidt coefficients c p , it is justified to keep 
only the D terms with coefficients of largest magnitude. Then, the contraction 
of the tensors A? and A® along the bond i—j with dimension D yields a good 
approximation to the true value M.if. 

[•^y'Lud r'u'd' ~ i^i] Ipud \Aj\ pr'u'd'' 

p=\ 

This method applied to all bonds with increased dimension provides us with the 
starting state for the optimization. 



6.4.1 Examples 

Let us now illustrate the variational methods with some examples. Models to 
which the PEPS algorithms have already been applied to include the Heisen- 
berg antiferromagnet [118] . the Shastry-Sutherland model [S2] and the system 
of hard-core bosons in a 2D optical lattice [75]. In the following, we recite the 
results for the latter system - which include calculations of ground state prop- 
erties and studies of the time-evolution after sudden changes in the parameters. 

The system of bosons in a 2D optical lattice is characterized by the Bose- 
Hubbard Hamiltonian 

H = - J (4 a i + h - c ) + J ni{hi - 1) + ^2 Vtfit, 

<i,j> i i 

where a\ and are the creation and annihilation operators on site i and 
hi = aja, is the number operator. This Hamiltonian describes the interplay 
between the kinetic energy due to the next-neighbor hopping with amplitude J 
and the repulsive on-site interaction U of the particles. The last term in the 



68 



-0.58 



-0.62 
E/JN 2 
-0.66 

-0.7 



0.7012 
0.7014 
0.7016 
0.7018 
-0.702 
0.7022 



exact 

- exact evolution 
■ Gutzwiller 

- D=2,3,4,5 




10 



10 



Figure 19: (taken from [75]) Energy as a function of time for the imaginary 
time-evolution of the system of hard-core bosons on a 4 x 4-lattice. The evolu- 
tions are performed sequentially with PEPS of virtual dimension D = 2, D = 3, 
D = 4 and D = 5. The times at which D is increased are indicated by verti- 
cal lines. For comparison, the exact ground state-energy, the exact imaginary 
time-evolution and the energy of the optimal Gutzwiller ansatz are included. 



Hamiltonian models the harmonic confinement of magnitude V{ = Vo(i — iq) 2 . 
Since the total number of particles N = J^, hi is a symmetry of the Hamiltonian, 
the ground-state will have a fixed number of particles. This number can be cho- 
sen by appending the term —\xN to the Hamiltonian and tuning the chemical 
potential fi. In the limit of hard-core interaction, U/J — > oo, two particles are 
prevented from occupying a single site. This limit is especially interesting in one 
dimension where the particles form the so-called Tonks-Girardeau gas [50] [HZ] ■ 
The particles in this gas are strongly correlated - which leads to algebraically 
decaying correlation functions. In two dimensions, the model was studied in 
detail in [58]. In the hard-core limit, the Bose-Hubbard model is equivalent to 
a spin-system with XX-interactions described by the Hamiltonian 

<i,j> i 
(i) (i) (i) 

Here, a x \ o\ and a z denote the Pauli-opcrators acting on site i. This Hamil- 
tonian has the structure that can be simulated with the PEPS algorithm: it 
describes L 2 physical systems of dimension d = 2 on a L x L-square lattice. 

In fig. [T5] the energy in the case of a 4 x 4-lattice is plotted as the system 
undergoes an imaginary time-evolution. Thereby, a time-step St = — i0.03 
is assumed and the magnitude of the harmonic confinement (in units of the 
tunneling-constant) is chosen as Vq/J = 36. In addition, the chemical potential 



69 




-0.31 



2 4 t 6 8 10 

Figure 20: (taken from [75]) Energy as a function of time for the imaginary 
time-evolution of the system of hard-core bosons on a 11 x 11-latticc. The 
evolutions are performed sequentially with PEPS of virtual dimension D = 2, 
D = 3, D = 4 and D = 5. The times at which D is increased are indicated by 
vertical lines. For comparison, the energy of the optimal Gutzwiller ansatz is 
included. 

is tuned to fi/J = 3.4 such that the ground state has particle-number (N) = 4. 
With this configuration, the imaginary time-evolution is performed both exactly 
and variationally with PEPS. As a starting state a product state is used that 
represents a Mott-like distribution with 4 particles arranged in the center of the 
trap and none elsewhere. The variational calculation is performed with D = 2 
first until convergence is reached; then, evolutions with D = 3, D = 4 and D = 5 
follow. At the end, a state is obtained that is very close to the state obtained by 
exact evolution. The difference in energy is \Eo=b — E exact \ ^ 6.4614 • 10 _5 J. 
For comparison, also the exact ground-state energy obtained by an eigenvalue- 
calculation and the energy of the optimal Gutzwiller ansatz are included in 
fig. HU The difference between the exact result and the results of the imaginary 
time-evolution is due to the Trotter-error and is of order 0(8t 2 ). The energy 
of the optimal Gutzwiller- Ansatz is well seperated from the exact ground-state 
energy and the results of the imaginary time-evolution. 

In fig. 1201 the energy as a function of time is plotted for the imaginary time- 
evolution on the 11 x 11-lattice. Again, a time-step St = — i0.03 is assumed 
for the evolution. The other parameters are set as follows: the ratio between 
harmonic confinement and the tunneling constant is chosen as Vq/J = 100 and 
the chemical potential is tuned to fj,/J = 3.8 such that the total number of 
particles (N) is 14. The starting state for the imaginary time-evolution is, 
similar to before, a Mott-like distribution with 14 particles arranged in the 
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Figure 21: (taken from [75]) Time evolution of the condensate density starting 
from a Mott-distribution with 14-particles arranged in the center of the trap. 
The magnitude of the trapping potential is Vq/J = 100. For the evolution, the 
Gutzwillcr ansatz and PEPS with D = 2, D = 3 and D = 4 are used. The inset 
shows the overlap between the D = 2 and D = 3-PEPS (solid line) and the 
D = 3 and D = 4-PEPS (dashed line). 

center of the trap. This state is evolved within the subset of PEPS with D = 2, 
D = 3, D = 4 and D = 5. As can be gathered from the plot, this evolution 
shows a definite convergence. In addition, the energy of the final PEPS lies well 
below the energy of the optimal Gutzwillcr ansatz. 

An application of the time-evolution algorithm with PEPS is found in the 
study of dynamic properties of hard-core bosons on a lattice of size 11 x 11. 
Here, the responses of this system to sudden changes in the parameters are 
investigated and the numerical results are compared to the results obtained by 
a Gutzwiller ansatz. An interesting property that is observed is the fraction of 
particles that are condensed. For interacting and finite systems, this property is 
measured best by the condensate density p which is defined as largest eigenvalue 
of the correlation-matrix (a\a,j). 

In fig. [2U the evolution of a Mott-distribution with 14 particles arranged in 
the center of the trap is studied. It is assumed that Vq/J — 100, fx/J = 3.8 and 
6t = 0.03. To assure that the results are accurate, the following procedure was 
used for simulating the time-evolution: first, the simulation has been performed 
using PEPS with D — 2 and D — 3 until the overlap between these two states 
fell below a certain value. Then, the simulation has been continued using PEPS 
with D = 3 and D = 4 as long as the overlap between these two states was 
close to 1. The results of this calculation can be gathered from fig. [21] What 
can be observed is that there is a definite increase in the condensate fraction. 
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Figure 22: (taken from [75]) Distance K between the time-evolved state and 
the state with reduced virtual dimension. The virtual dimensions D = 2, D = 3 
and D = 4 arc included. The distance is plotted for the evolution of a Mott- 
distribution with ./V = 14, as explained in fig. [21] From the inset, the deviation 
of the particle number from the value 14 can be gathered. 

The Gutzwillcr ansatz is in contrast to this result since it predicts that the 
condensate density remains constant. The inset in fig. [5U shows the overlap of 
the D = 2 with the D = 3-PEPS and the D = 3 with the D = 4-PEPS. 

Finally, we make a few comments about the accuracy of the algorithm. One 
indicator for the accuracy is the distance between the time-evolved state and 
the state with reduced virtual dimension. For the time-evolution of the Mott- 
distribution that was discussed before, this quantity is plotted in fig. [22] We find 
that the distance is typically of order 10 -3 for D = 2 and of order 10 -4 for D = 3 
and D = 4. Another quantity that is monitored is the total number of particles 
(N) . Since this quantity is supposed to be conserved during the whole evolution, 
its fluctuations indicate the reliability of the algorithm. From the inset in fig. l22l 
the fluctuations of the particle number in case of the time-evolution of the 
Mott-distribution can be gathered. We find that these fluctuations are at most 
of order 10 -5 . 

6.4.2 PEPS and fermions 

The critical reader should by now have complained that we are only talking 
about spin systems but not about fermionic systems. Indeed, one of the long 
term goals of the numerical approaches discussed here is to be able to simulate 
e.g. the Fermi-Hubbard model in the relevant parameter regime. 

The methods that we discussed in 1 dimension are perfectly applicable to 
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fcrmionic systems, as the amazing Jordan- Wigner transformation allows to map 
a local Hamiltonian of fcrmions to a local Hamiltonian of spins, and we know 
that the MPS-techniques work provably well on the latter. The big problem 
however is that the Jordan- Wigner transformation only works in 1 dimension: 
if we use it on a 2-D lattice, a local Hamiltonian of fcrmions is mapped to a 
highly nonlocal Hamiltonian of spins. PEPS on the other hand are devised to be 
such that they have extremal local properties; if the Hamiltonian contains a lot 
of strong nonlocal terms, we can not expect a PEPS to exhibit the correspond- 
ing extremal long-range correlations. The natural question to ask is therefore 
whether there exists a generalization of the Jordan- Wigner transformation to 
higher dimensions. This was indeed shown to be possible in |115) : given any 
local Hamiltonian in terms of fermionic operators such as the Hubbard model 
in 2 or 3 dimensions, then there exists a local spin 3/2 Hamiltonian whose low- 
energy sector corresponds exactly to the original fcrmionic Hamiltonian. The 
conclusion is that the PEPS methods are equally applicable to quantum spin 
systems as to fcrmionic systems. 

Another and more efficient approach is to make use of quantum numbers. 
In general, it is difficult to keep track of quantum numbers on a 2-D lattice. 
An exception however is given by the parity of the occupation number between 
two sites: by blocking the PEPS tensors in a specific way, one can invoke the 
even or odd occupation number between any sites and as such eliminate the fact 
that sprurious effective long-range interactions arise in fermionic lattice systems 
[3"0] . This is very relevant as it leads to a huge speed-up of the algorithms for 
fermionic lattice systems. 

6.5 PEPS on infinite lattices 

In parallel with the 1-dimensional MPS, we can also explicitly make use of the 
translational symmetry in the class of PEPS to simulate low-energy properties 
of infinite lattices. A naive procedure is to impose complete translational sym- 
metry and do a line-search (such as conjugate gradient) within the parameter 
space of the tensor; note that this can be done as expectation values of the corre- 
sponding PEPS can be calculated using the appropriate 1-D infinite algorithms 
discussed above. However, the cost function is highly nonlinear and the number 
of parameters grows very fast with D, such that those brute-force methods will 
typically lead to local minima. A smarter approach is to use imaginary time 
evolution, but keeping the translational invariance. It was already shown how 
to do this in the 1-D case, where a particular kind of Trotter expansion lead to a 
translational invariant matrix product operator, and we can repeat this in 2-D 
leading to a translational invariant PEPS-operator. The only nontrivial part 
is the question of how to reduce the bond dimension after one time evolution 
step; but here we can again get inspiration of the 1-D case, and ask the question 
which projectors we can put on the bonds such as to maximize the overlap of 
the projected one with the evolved one. 

A more straightforward approach is to assume an ..ABAB... symmetry where 
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Figure 23: (taken from [56]) Magnetization m x (h) in the ground state \ ^h) 
of the two-dimensional quantum Ising model with transverse magnetic field. A 
trivial iPEPS (inner dimension D — 1) produces a magnetization that grows es- 
sentially linearly with the magnetic field until it saturates. Instead, the simplest 
non-trivial iPEPS (inner dimension D = 2) produces a magnetization curve that 
overlaps with the results of series expansions for both small and large magnetic 
fields. [D = 3 leads to results that could hardly be distinguished from those 
for D = 2 in this figure]. Notice that around h « 3.1 the derivative of the 
magnetization m x (h) changes suddenly. 
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each tensor A has all its neigbouring tensors B and vice-versa (of course this 
is only possible for bipartite lattices; a different choice can be made for differ- 
ent types of lattices). Again, imaginary time evolution can be used to find the 
ground state, and this was first studied in the paper [56] where the term iPEPS 
(infinite PEPS) was coined. The idea is as follows: take the even-odd-horizontal- 
vertical Trotter decomposition as discussed previously, and next evolve with just 
one operator acting on two nearest neigbour sites. Effectively, this increases 
the bond dimension between those two sites (A and B). The environment of 
those spins can be readily calculated using the infinite 1-D translational invari- 
ant methods discussed above, and then the variational problem becomes the 
problem of finding new A' and B' that approximate the one with higher bond 
optimally. Again, this can be done using the alternating least squares method. 
Subsequently, we replace all tensors A and B with A 1 and B' , and continue until 
convergence. The last step - replacing all tensors with the optimal local ones - 
is only justified if the time step in the imaginary time evolution is very small, 
but in practice, this seems to work very well, as illustrated in figure [531 
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7 Conclusion 



Recent progress in Quantum Information Theory has provided scientists with 
new mathematical tools to describe and analyze many-body quantum systems. 
These new tools have given rise to novel descriptions of the quantum states 
that appear in Nature, which are very efficient both in terms of the number of 
variables used to parametrize states and the number of operations to determine 
expectation values of typical observables. In a sense, they allow us to describe 
the "corner of Hilbert space" where relevant states are located with an effort 
that only scales polynomially with the number of particles, as opposed to the 
exponential scaling resulting with other descriptions. These results have auto- 
matically led to a diverse set of new powerful algorithms to simulate quantum 
systems. Those algorithms allow us to describe ground states, thermal equilib- 
rium, low excitations, dynamics, random systems, etc, of many-body quantum 
systems, and thus to attack new kinds of problems obtaining very precise re- 
sults. Moreover, the methods work in one and more spatial dimensions. In 
the first case, the success of some of those methods is directly related to the 
extraordinary performance of DMRG. In higher dimensions, they also give rise 
to a better understanding of several many-body systems for which a description 
has not been possible with the existing techniques. 

This paper has reviewed these new methods in a unified form. We have 
introduced MPS and its extensions to higher dimensions, PEPS, and shown 
how one can build powerful algorithms that find the best descriptions of states 
within that family of states. The algorithms are relatively simple to implement, 
although they require some tricks that have been reported in this paper as well. 
We have also given simple matlab codes in Appcndix[C]to illustrate how one can 
program some of the methods. Thus, we believe that the present paper may be 
very useful both to scientist interested in implementing these new algorithms to 
describe any kind of many-body quantum systems, as well as those interested 
in creating new algorithms for some specific purposes. We have also provided 
several appealing evidences of the fact that most interesting states in Nature are 
well described by MPS and PEPS. This, in fact, indicates that our algorithms 
as well as future extensions can provide us with unique tools to explore the 
fascinating physics of quantum many-body systems. 
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A Local reduced density operators of spin sys- 
tems 

Due to the variational nature of ground states, there always exists a ground state 
with the same symmetries as the associated Hamiltonian. If the Hamiltonian 
has translational symmetry and consists of 2-body nearest neighbor interac- 
tions, then it is clear that the energy of a state with the right symmetry is 
completely determined by its reduced density operator of 2 neighboring spins. 
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Figure 24: (taken from [116] ) Convex sets of the possible reduced density op- 
erators of translational invariant spin 1/2 states in the XX-ZZ plane. The big 
triangle represents all positive density operators; the inner parallellogram rep- 
resents the separable states; the union of the separable cone and the convex hull 
of the full curved line is the complete convex set in the case of a 1-D geometry, 
and the dashed lines represent extreme points in the 2-D case of a square lattice. 
The singlet corresponds to the point with coordinates (—1,-1). 

The reduced density operators arising from these (eventually mixed) states with 
a given symmetry form a convex set, and the energy for a given Hamiltonian will 
be minimized for a state whose reduced density operator is an extreme point 
in this set. More specifically, the equation Tr(Hp) = E determines a hyper- 
plane in the space of reduced density operators of states with a given symmetry, 
and the energy will be extremal when the hyperplane is tangent to the convex 
set (E = E extr ). The problem of finding the ground state energy of nearest 
neighbor translational invariant Hamiltonians is therefore equivalent to the de- 
termination of the convex set of 2-body reduced density operators arising from 
states with the right symmetry. Strictly speaking, these two problems are dual 
to each other. In the case of quadratic Hamiltonians involving continuous vari- 
ables, the determination of this convex set was solved for fairly general settings 
in [143] by means of Gaussian states. The determination of this convex set in 
the case of spin systems however turns out to be much more challenging. 
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Let us illustrate this with a simple example. 

H = - s ? S j + S i S j + AS ! S j 

<i,j> 

on a lattice of arbitrary geometry and dimensional. Due to the symmetries, the 
reduced density operator of two nearest neighbors can be parameterized by only 
two parameters l 27 l: 

P = - (/ g) / + X(a x ® d X + dy ® <j y ) + z<r z ® a z ) . 

Positivity of p enforces — 1 < z <1 — 2|x|, and the state is separable iff 1 + z < 
2\x\. In the case of an infinite 1-D spin chain, the ground state energy E(A) 
has been calculated exactly [18] . and this determines the tangent hyperplanes 

2x + zA + E(A) = 

whose envelope makes up the extreme points of the convex set of reduced density 
operators of translationally invariant 1-D states: the boundary of this convex 
set is parameterized by 

z = -dE(A)/dA 

x = -(E{A)+dE(A)/8A)/2, 

which we plotted in Figure [M] We also plot the boundary for the 2-dimensional 
square lattice. These 2-D data were obtained by numerical methods [1 181 1 1 1 11 
[70]); of course this convex set is contained in the previous one, as all the semidcf- 
inite constraints defining the set corresponding to 1-D are strictly included in 
the set of constraints for the 2-D case. Finally, we plot the set of separable 
states, which contains the reduced density operators of the allowed states for 
a lattice with infinite coordination number. The boundary of this separable 
set is given by the inner diamond; this immediately implies that the difference 
between the exact energy and the one obtained by mean field theory will be 
maximized whenever the hyperplane that forms the boundary of the first set 
will be parallel to this line. This happens when A = — 1 (independent of the 
dimension!), which corresponds to the antiferromagnetic case, and this proves 
that the " entanglement gap" [27] in the XXZ-plane is maximized by the antifer- 
romagnetic ground state for any dimension and geometry. Similarly, it proves 
that the ground state is separable whenever A > 1 and A = — oo. Note also 
that in the 2-D case, part of the boundary of the convex set consists of a plane 
parameterized by 2x + z + E(l) = 0. This indicates a degeneracy of the ground 
state around the antiferromagnetic point, and indicates that a phase transition 
is occurring at that point (more specifically between an Ising and a Berezinskii- 
Kostcrlitz-Thouless phase) . 

26 We assume that the graph corresponding to the lattice is edge-transitive, meaning that 
any vertex can be mapped to any other vertex by application of the symmetry group of the 
graph. 

27 This can easily be proven by invoking local twirling operations which leave the Hamilto- 
nian invariant. 
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Figure 25: (taken from [116] ) Convex sets in the XXZ-plane: the inner diamond 
borders the set of separable states (see Fig. 1). Dash-dotted: extreme points of 
the convex set produced by MPS of D = 2. 

As a good illustration of the actual accuracy obtained with MPS, we cal- 
culated the convex set obtained with MPS in the thermodynamic limit for the 
XXZ-chain with D = 2, where D is the dimension of the matrices in the MPS 
(see figure l2"5j) . It is almost unbelievable how good the exact convex set can be 
approximated. Note that typical DMRG calculations have D ~ 200, and that 
the accuracy grows superpolynomial in D. Note also that the D = 1 case cor- 
responds to mean-field theory, whose corresponding convex set coincides with 
the set of separable states. 

The same argument involving the notion of a correlation length applies in 
higher dimensions and indicates that PEPS represent ground states of gapped 
local Hamiltonians well. Note however that the convex set in the 2-D case 
is much closer to the separable one than in the 1-D case; this gives a hint 
that PEPS of smaller dimension will suffice to give the same accuracy as in 
the 1-D case. In the next section we will quantitatively bound how well a 
translationally invariant state can be represented in terms of a MPS, and will 
analyze the corresponding implications for the description of ground states of 
ID spin chains. 
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B MPS represent ground states faithfully 

We will derive an upper bound to the error made by approximating a general 
ground state of a 1-D quantum spin system by a MPS. As we will show below, 
this has very important implications in the performance of the renormalization 
algorithms to describe ground states of ID spin chains. 

Lemma 1 There exists a MPS \ipE>) of dimension D such that 

N-i 

|||V> -rfe}|| 2 < 2 $> a p) 

c*=l 

where e a (D)=Y."= D+ i V [a]l - 

Proof: We can always write ip as a MPS of dimension D = 2 N ' 2 and fulfilling 
V^^Mi^M't = j ^4[™]«t 7 Y[" i + 1 ] J 4[" l ] i — j\l m ] 

i i 

Let us now consider the D-dimensional MPS \4>d) which is defined by the Dx D 
matrices [^4^ q '*](i...d.i...l>) (i-e. the upper-left block oi A^ 1 ). The goal is now 
to bound (tpltpu). The gauge conditions were chosen such as to make the task 
simple: 

(M4>) =Tr[$ 2 (...$ Ar _ 2 ($ JV _ 1 (A Ar - 1 P)P)P...)P] ; (32) 

here P = Y.k=i \ k ) ( k \ and $m(X) = Ei 4 Mit 14 Wi represents a trace-preserving 
completely positive map (TPCP-map) parameterized by the Kraus operators 
^4. ["*]*. Let us now recursively define 

y[k] = $fc (y[k + l] p ^ ; y[N-l] = A [N-l]p. 

observe that A^ = $ fe (A[ fc+1 l). We want a bound on Tr|AM — Yi|, as equation 
(|32|) is equal to Tr(Y^). The crucial property wc need is that TPCP-maps are 
contractive with relation to the trace-norm! 28 !: Tr|$(X)| < Tr|X|. It follows that 

Tr|AM - yH| = Tr |$ k (A^ - y^p) | < 

< Tr | A [k+i] _ Y [ k+1 lp| 

< Tr|A[ k+1 l - Y[ k+1 l| + Tr|Al k+1 ](I - P)|. 

ON/2 r, , 

Note that the last term in the sum is exactly given by E q =d+i ^ ■ The 
theorem now follows immediately by recursion and by observing that (iI>d\iPd) < 
1 by similar arguments. □ 

28 This can directly be proven by considering the Neumark representation of a TPCP-map 
as a unitary in a bigger space. 
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The implications of this result are very strong: it shows that for systems 
for which the e a (D) decay fast in D, there exist MPS with small D which 
will not only reproduce well the local correlations (such as energy) but also 
all the nonlocal properties (such as correlation length). The following lemma 
now relates the derived bound to the Renyi entropies of the reduced density 
operators, through which one can make the connection to the ground states of 
ID Hamiltonians. The Renyi entropies of p arc defined as 

S a (p) = -^—log(Trp a ), 
1 — a 

and we will consider < a < 1. We denote as before e(D) = Y^=D+i ^» with 
Xi the nonincreasingly ordered eigenvalues of p. Then we have 

Lemma 2 Given a density operator p. If < a < 1, then log(e(D)) < 
(<?«(/>)- log 

Proof: Let us first characterize the probability distribution that has maximal 
possible weight in its tail (i.e. p = J2i^=D+iPi) f° r a gi yen Renyi-entropy. 
Introducing a free parameter < h < (l—p)/D, such a probability distribution 
must be of the form 

Pi = l-p-(D-l)h 

h = P2 = P3 = ■ ■ -PD+p/h 
PD+p/h+li ■ ■ 'Poo = 

because this distribution majorizes all other ones with given p,D,pjj (Rcnyi- 
cntropies are Schur-convex functions). For a given p, D, h, it holds that 

Y,Pi - (l-p-(D-l)h) a + (D-l+p/h)h a 

i 

> D^+ph - 1 . 
Minimizing this expression with relation to h, we get 

^p^rVl/ffl-a) 1 ^). 

i 

Denoting S a (p, D) the minimal possible entropy for given p, D, we get 

1 / D^ a p a 

S a (p,D)>- log' 



1-Q' V(l-a) 1 -«a 

and hence 

P < cxp t 1 —^- (s a (p, D) - log -5- 
\ a \ 1 — a, 

The proof now follows by replacing S a (p, D) by S a (p) . □ 
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This lemma is very interesting in the light of the fact that in the case of 
critical systems, arguable the hardest ones to simulate Q the Renyi-entropy of 
a contiguous block of L spins scales as [BH [S3 HH 1131] 



for all a > 0; here c is the central charge. The fact that the eigenvalues of 
Pl decay fast has previously been identified as a indication for the validity of 
the DMRG-approach [89] . The truncation error |137( [97] , which has been used 
in the DMRG community as a check for convergence, is essentially given by 
e(D) — e(2D) and therefore indeed gives a good idea of the error in a simulation. 

Let us investigate how the computational effort to simulate such critical 
systems scales as a function of the length N = 2L of the chain. Let us therefore 
consider the Hamiltonian associated to a critical system, but restrict it to 2L 
sites. The entropy of a half chain (we consider the ground state \i[>ex) of the 
finite system) will typically scale as in eq. (|33| but with an extra term that 
scales like 1/7V. Suppose we want to enforce that HlV'ex) — l^rOH 2 < £o/L with 
eo independent of L £3 Denote the minimal D needed to get this precision for 
a chain of length 2L by Dl,. Following lemma (1) and the fact that the entropy 
of all possible contiguous blocks reaches its maximum in the middle of the chain 
(hence p < eo/L 2 is certainly sufficient), lemma (1) and (2) combined yield 



This shows that D only has to scale polynomially in L to keep the accuracy 
cq/L fixed; in other words, there exists an efficient scalable representation for 
ground states of critical systems (and hence also of noncritical systems) in terms 
of MPS! Such a strong result could not have been anticipated from just doing 
simulations. Furthermore, Hastings has proven that ground states of gapped 
systems always obey a strict area law [35] ■ This implies that the ground state 
of any gapped spin chain is indeed well approximated by a MPS. For a more 
detailed description of the relations between area laws and approximability, we 
refer to [TOO] . 

Now what about the complexity of finding this optimal MPS? It has been 
observed that DMRG converges exponentially fast to the ground state with a 
relaxation time proportional to the inverse of the gap A of the system [97] . 
For translational invariant critical systems, this gap seems to close only poly- 
nomially. As we have proven that D only have to scale polynomially too, the 
complexity of deciding whether the ground state energy of TD quantum systems 

29 For non— critical systems, the renormalization group flow is expected to increase the Renyi 
entropies in the UV direction. The corresponding fixed point corresponds to a critical system 
whose entropy thus upper bounds that of the non-critical one. 

30 We choose the 1/L dependence such as to assure that the absolute error in extensive 
observables does not grow. 




(33) 
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is below a certain value is certainly in NP (as it can be checked efficiently if Mer- 
lin gives the MPS-description). The problem would even be in P if the following 
conditions are met: 1) the a-entropy of blocks in the exact ground state grow 
at most logarithmically with the size of the block for some a < 1; 2) the gap of 
the system scales at most polynomially with the system size; 3) given a gap that 
obeys condition 2, there exists an efficient DMRG-likc algorithm that converges 
to the global minimum. As the variational MPS approach [121) is essentially an 
alternating least squares method of solving a non-convex problem which is in 
worst case NP-hard [55] , there is a priori no guarantee that it will converge 
to the global optimum, although the occurrence of local minima seems to be 
unlikely [97| . Surprisingly, one can indeed construct a family of Hamiltonians 
with nearest neighbour interactions on a line for which the ground state is an 
exact MPS with polynomial D, but for which it is an NP-complete problem to 
find it [99] ■ As the corresponding Hamiltonian has a gap that closes polynomi- 
ally in the size of the system, the only hope that VMPS/DMRG methods to be 
in P is in the case of gapped systems 0. 

C Mat lab code 

As a last part, we would like to give an idea of how to program the variational 
methods explained in the previous sections. We present two functions, one for 
the calculation of ground- and first excited states and one for the reduction of 
the virtual dimension of matrix product states. We demonstrate these functions 
by means of the antifcrromagnctic Hcisenberg chain. 

C.l Minimization of the Energy 

The function minimizeE optimizes the parameters of a matrix product state 
in such a way that the expectation value with respect to a given Hamiltonian 
tends to a minimum. The function expects this Hamiltonian to be defined in a 
M x N cell hset, where N denotes the number of sites and M the number of 
terms in the Hamiltonian. Assuming the Hamiltonian is of the form 

M 
m— 1 

31 Let us specify an alternative method which should in principle not get trapped in local 
minima in the case of gapped systems. Like in the adiabatic theorem, we can construct a 
time dependent Hamiltonian H(t) with H(0) trivial and -ff(l) the Hamiltonian to simulate; 
if we discretize this evolution in a number of steps that grows polynomially in the inverse 
gap, the adiabatic theorem guarantees that we will end up in the ground state of H(l) if we 
can follow the ground state of H(t) closely. The idea is to make D of \ip£>(t)) large enough 
such as to follow the ground state |^>(t)) close enough in such a way that the optimization 
is always convex around the global optimum within the domain — < e. As we 

are simulating this classically, we could even do imaginary time evolution over a longer time 
during each step. 
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the element hset{m , j } equals Km . Further arguments are the virtual dimension 
of the resulting matrix product state, D, and the expected accuracy of the energy, 
precision. 

Output arguments are the optimized energy E and corresponding matrix 
product state mps. The matrix product state is stored as a 1 x TV cell, each 
entry corresponding to one matrix. 

Optionally, a matrix product state mpsB can be specified as an argument 
to which the resulting state shall be orthogonal. This is especially useful for 
calculating the first excited state. 

function [E, mps] =minimizeE(hset,D, precision, mpsB) 

[M,N]=size(hset) ; 
d=size(hset-(l,l},l); 
mps=cr eater andommps (N , D , d) ; 
mps=prepare (mps) ; 

7, storage-initialization 
Hstorage=initHstorage (mps , hset , d) ; 

if "isempty (mpsB) , Cstorage=initCstorage (mps , [] ,mpsB , N) ; end 
P= []; 

7. optimization sweeps 
Evalues= [] ; 

■/, ****************** cycle 1: j -> j+1 (from 1 to N-l) **************** 
for j=i: (N-l) 

7i projector-calculation 
if "isempty (mpsB) 
B=mpsB{j>; 
Clef t=Cstorage-t j } ; 
Cright=Cstorage{j+l>; 

P=calcprojector_onesite(B J Cleft J Cright) ; 

'/, optimization 
Hleft=Hstorage(: ,j) ; 
Hright=Hstorage ( : , j+1) ; 
hsetj=hset(: ,j) ; 

[A, E] =minimizeE_onesit e (hset j , Hlef t , Hright , P) ; 
[A,U]=prepare_onesite(A, J lr') ; 
mps{j}=A; 

Evalues= [Evalues,E] ; 

'/, storage-update 
for m=l:M 

h=reshape (hset-tm , j } , [1 , 1 , d , d] ) ; 

Hstorage{m,j+l>=updateCleft(Hleft-[m> J A J h J A) ; 

if "isempty (mpsB) 

Cstorage-tj+l}=updateCleft(Cleft J A J [] ,B) ; 

% ****************** cycle 2: j -> j-1 (from N to 2) ****************** 
for j=N: (-1) :2 

'/, projector-calculation 
if "isempty (mpsB) 
B=mpsB{j>; 
Clef t=Cstorage-t j } ; 
Cright=Cstorage{j+l>; 

P=calcprojector_onesite (B , Cleft , Cright ) ; 

Hlef t=Hstorage ( : , j ) ; 
Hright=Hstorage ( : , j+1) ; 
hsetj=hset(: ,j) ; 

[A , E] =minimizeE_onesit e (hset j , Hlef t , Hright , P) ; 
[A,U]=prepare_onesite(A, 'rl') ; 
mps{j}=A; 

Evalues= [Evalues,E] ; 

7 t storage -update 
for m=i:M 

h=reshape(hset{m J j>, [l.l.d.d]) ; 

Hstorage-tm, j }=updateCright (Hright{m>, A,h,A) ; 

if "isempty (mpsB) 

Cstorage-tj}=updateCright (Cright , A, [] ,B) ; 
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if (std(Evalues) /abs (mean(Evalues) ) precision) 
mps{l}=contracttensors Cmps-tl} , 3 , 2 , U, 2 , 1) ; 
mps-[l>=perinute(mps-[l}, [1,3,2] ) ; 

function [A , E] =minimizeE_onesite (hset j , Hlef t , Hright , P) 

DAl=size (Hlef t{l> , 1) ; 
DAr=size (Hright{l> , 1) ; 
d=size(hsetj{l},l) ; 

'/, calculation of Heff 
M=size(hsetj ,1) ; 

Hef f =0 ; 
for m=l:M 

Hef f m=contracttensors (Hlef t{m> ,3,2, Hright {m} ,3,2); 
Heffm=contracttensors(Heffm,5,5,hsetj-[in},3,3) ; 
Heffm=permute(Heffm, [1 ,3 ,5 ,2 ,4,6] ) ; 
Heffm=reshape(Heffm, [DAl*DAr*d,DAl*DAr*d] ) ; 
Heff=Heff+Heffm; 

'/, projection on orthogonal subspace 
if ~isempty(P), Hef f =P ' *Hef f *P ; end 

7, optimization 
options . disp=0 ; 

[A,E] =eigs (Heff , 1, »«r* .options) ; 
if 'isempty(P) , A=P*A; end 
A=r eshape (A, [DAI , DAr , d] ) ; 



function [P] =calcprojector_onesite (B , Cleft ,Cright) 

y=contracttensors (Clef t ,3 ,3 ,B ,3 , 1) ; 
y=contracttensors(y,4, [2,3] ,Cright,3, [2,3] ) ; 
y=permute(y, [1,3,2]) ; 
y=r eshape (y, [prod ( size (y) ) ,1] ) ; 

Q=orth([y,eye(size(y,l))]); 
P=q(: ,2:end) ; 



C.2 Time Evolution 

The function reduceD forms the basis for the simulation of a time evolution. It 
multiplies a given matrix product state with a given matrix product operator 
and reduces the virtual dimension of the rcsuling state, i.e. it searches a ma- 
trix product state with reduced virtual dimension and minimal distance to the 
original state. The matrix product state and the matrix product operator are 
specified in the arguments mpsA and mpoX. As before, they are represented by 
a cell with entries identifying the matrices. The reduced virtual dimension is 
specified in the argument DB. The argument precision defines the convergence 
condition: if fluctutions in the distance are less than precision, the optimiza- 
tion is assumed to be finished. 

The output argument is the optimized matrix product state mpsB with vir- 
tual dimension DB. 

function mpsB=reduceD (mpsA,mpoX,DB , precision) 

N=length(mpsA) ; 

d=size(mpsA-[l},3) ; 

mpsB=cr eater andommps (N , DB , d) ; 

mpsB=prepare (mpsB) ; 

'/, initialization of the storage 

Cstorage=initCstorage(mpsB,mpoX,mpsA,N) ; 
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7, optimization sweeps 
Kvalues= [] ; 

■/, ****************** cycle 1: j -> j+1 (from 1 to N-l) **************** 
for j=i: (N-l) 

7. optimization 

Clef t=Cstorage{j } ; 

Cright=Cstorage{j+l>; 

A=mpsA-tj>; X=mpoX{j>; 

[B,K]=reduceD2_onesite (A, X, Cleft, Cright) ; 
[B,U]=prepare_onesite(B, 'lr') ; 
mpsB{j}=B; 

Kvalues= [Kvalues,K] ; 
7. storage-update 

Cstorage{j+l}=updateClef t(Clef t,B,X,A) ; 

•/, ****************** cycle 2: j -> j-1 (from N to 2) ****************** 
for j=N: (-1) :2 

7. optimization 

Clef t=Cstorage{ j } ; 

Cright=Cstorage{j+l>; 

A=mpsA{j}; X=mpoX{j>; 

[B J K]=reduceD2_onesite(A,X,Cleft,Cright) ; 
[B,U]=prepare_onesite(B, 'rl') ; 
mpsB-tj>=B; 

Kvalues= [Kvalues.K] ; 
7. storage-update 

Cstorage{j}=updateCright(Cright,B,X,A) ; 

if std (Rvalues) /abs (mean (Rvalues) ) <precision 
mpsB{l}=contracttensors(mpsB{l},3,2,U,2,l) ; 
mpsB{l}=permute(mpsB{l}, [1,3 ,2] ) ; 



function [B , K] =reduceD2_onesite (A , X , Clef t , Cright ) 

Clef t=contracttensors (Clef t, 3 ,3, A, 3,1} ; 

Clef t=contracttensors (Cleft, 4, [2 ,4] ,X,4, [1 ,4] ) ; 

B=contracttensors(Cleft,4, [3,2] , Cright, 3, [2,3]) ; 
B=permute(B, [1,3,2]); 

b=reshape(B, [prod (size (B) ) , 1] ) ; 
K=-b J *b; 



C.3 Auxiliary functions 

The previous two functions depend on several auxiliary functions that are printed 
in this section. 

• Gauge transformation that prepares the MPS mpsB in such a form that 
N e ff is equal to the identity for the first spin (see section [331) : 

function [mps] =prepare (mps) 
N=length(mps) ; 
for i=N:-l:2 

[mps{i>,U]=prepare_onesite(mps{i>, 'rl') ; 

mps{i-l}=contracttensors (mps-ti-1} ,3 ,2 ,U,2 , 1) ; 
mps{i-l}=permute(mps-(i-l}, [1,3,2] ) ; 



function [B,U,DB]=prepare_onesite( A, direction) 

[Dl,D2,d]=size(A) ; 
switch direction 
case 'lr J 

A=permute(A, [3,1,2]) ; A=reshape(A, [d*Dl,D2] ) ; 
[B,S,U]=svd2(A) ; DB=size(S , 1) ; 
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B=reshape(B, [d.Dl.DB]) ; B=permute(B, [2,3,1]) ; 
U=S*U; 
case 'rl' 

A=permute(A, [1,3,2]) ; A=reshape(A, [Dl,d*D2] ) ; 
[U,S,B]=svd2(A) ; DB=size (S , 1) ; 
B=reshape(B, [DB,d,D2] ) ; B=permute(B, [1,3,2]); 
U=U*S ; 



• Initialization of storages: 



function [Hstorage] =initHstorage Cmps , hset , d) 

[K,N] =size(hset); 
Hstorage=cell(M,N+l) ; 

for m=l:H, Hstorage{m, 1>=1; Hstorage{m,N+l}=l; end 
for j=N:-l:2 
for m=i:M 

h=reshape(hset-[m J j> l [l.l.d.d]); 

Hstorage-(m, j }=updateCright (Hstorage{m, j + 1} ,mps-(j } ,h,mps-[j }) ; 



function [Cstorage] =initCstorage (mpsB ,mpoX,mpsA,N) 

Cstorage=cell(l,N+l) ; 

Cstorage{l}=l ; 

Cstorage{N+l}=l; 

for i=H:-l:2 

if isempty (mpoX) , X= [] ; else X=mpoX{i}- ; end 
Cstorage-ti}=updateCright (Cstorage{i+l} ,mpsB-(i} , X,mpsA-[i}) ; 



function [Cleft] =updateCleft (Cleft ,B,X, A) 

if isemptyCX), X=reshape (eye (size (B , 3) ) , [1,1,2,2]); end 

Clef t=contracttensors (A , 3 , 1 , Cleft ,3,3); 
Cleft=contracttensors(X,4, [1,4] , Cleft, 4, [4,2]) ; 
Cleft=contracttensors(conj(B) ,3, [1,3] , Cleft, 4, [4,2] ); 



function [Cright] =updateCright (Cright, B,X, A) 

if isemptyCX), X=reshape (eye (size (B ,3) ) , [1 , 1 ,2 ,2] ) ; end 

Cright=contracttensors (A , 3 , 2 , Cright ,3,3); 
Cright=contracttensors(X,4, [2,4] , Cright ,4, [4,2] ) ; 
Cright=contracttensors(conj (B) ,3, [2,3] .Cright ,4, [4,2] ) ; 



• Creation of a random MPS: 



function [mps] =creat erandommps (N , D , d) 
mps=cell(l,N) ; 

mps{l>=randn(l,D,d)/sqrt(D) ; 
mpsOJ>=randn(D , 1 ,d) /sqrt (D) ; 
for i=2: (SJ-1) 

mps{i>=randn(D,D,d)/sqrt(D) ; 

• Expectation value of the MPS mps with respect to the operator defined 
in hset: 

function [e ,n] =expectationvalue (mps , hset) 

[M,N]=size(hset) ; 
d=size(mps{l},3) ; 

'/, expectation value 

for m=l:M 

for j=N:-l:l 

h=hset{m, j>; 
h=reshape(h, [l.l.d.d]); 
em=updateCright(em,mps-[j},h,mps{j}) ; 
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11=1; 

X=eye(d) ; X=reshape(X, [1, l.d.d] ) ; 
for j=N:-l:l 

n=updateCright (n,mps-(j } ,X,mps-[j }) ; 

• Contraction of index indX of tensor X with index indY of tensor Y (X 
and Y have a number of indices corresponding to numindX and numindY 
respectively): 

function [X,numindX]=contracttensors(X,numindX,indX,Y, numindY, indY) 

Xsize=ones(l,numindX) ; Xsize(l : length(size(X) ) ) =size (X) ; 
Ysize=ones(l, numindY) ; Ysize(l : length(size(Y) ) )=size (Y) ; 

indXl=l:numindX; indXl(indX)=[] ; 
indYr=l:numindY; indYr(indY)= [] ; 

sizeXl=Xsize(indXl) ; 
sizeX=Xsize(indX) ; 
sizeYr=Ysize(indYr) ; 
sizeY=Ysize(indY) ; 

if prod(sizeX)~=prod(sizeY) 

error('indX and indY are not of same dimension.'); 

if isempty (indYr) 

if isempty (indXl) 

X=permute(X, [indX] ) ; 
X=reshape(X, [1 ,prod(sizeX)] ) ; 

Y=permute(Y, [indY]) ; 
Y=reshape(Y, [prod(sizeY) , 1] ) ; 

X=X*Y; 
Xsize=l; 



else 

X=permute(X, [indXl , indX] ) ; 

X=reshape(X, [prod(sizeXl) ,prod(sizeX)] ) ; 

Y=permute(Y, [indY]); 
Y=reshape(Y, [prod(sizeY) , 1] ) ; 

X=X*Y; 

Xsize=Xsize(indXl) ; 
X=reshape(X, [Xsize.l]) ; 

X=permute (X , [indXl , indX] ) ; 

X=reshape(X, [prod(sizeXl) .prod(sizeX) ] ) ; 

Y=permute (Y , [indY , indYr] ) ; 

Y=reshape(Y, [prod(sizeY) .prod(sizeYr)] ) ; 

X=X*Y; 

Xsize=[Xsize(indXl) .Ysize(indYr)] ; 
numindX=length(Xsize) ; 
X=reshape (X , [Xsize , 1] ) ; 



• Economical singular value decomposition: 



function [U J S J V]=svd2(T) 
[m,n] =size(T); 

if m>=n, [U,S,V]=svd(T,0) ; else [V, S,U] =svd(T ' , 0) ; end 
V=V ; 
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Figure 26: (Left) Error of the variational method as a function of the virtual 
dimension D for N ~ 10 spins. The blue triangles represent the ground state, 
the red triangles the first excited state. (Right) Time evolution as described 
in example 2. The magnetization of the central spin is represented by blue 
triangles and the magnetization of the spin adjacent to the central spin by red 
triangles. For comparison, exact results are also included (black lines). 



C.4 Examples 

As a first example, we show how to calculate the ground-state and the first ex- 
cited state of the antiferromagnetic Heisenberg chain using the method minimi zeE 
from before. In this example, the chain-length N is assumed to be 10 and the 
virtual dimension D is set to 5. 

K=10; 
D=5; 

precision=le-5 ; 

'/, Heisenberg Hamiltonian 
M=3*(N-1) ; 
hset=cell(M,N) ; 

BX=[0,i;l,0] ; By=[0,-li;ii,Q]; sz=[l ,0;0,-l] ; id=eye{2); 
for m=l : M, for j=l :N , hset{m , j}=id ; end; end 
for j=i! (H-l) 

hset{3*(j-l)+l,j}=sx; hset{3*(j-l)+l, j+l}=sx; 

hset{3*(j-l)+2,j}=sy; hset{3*(j-l)+2, j+l}=sy; 

hset{3*(j-l)+3,j}=sz; hset{3*(j-l)+3, j+l}=sz; 

% ground state energy 
randn( J state ' ,0) 

[E0 ,mps0] =minimizeE(hset ,D , precision, [] ) ; 
fprintfC'EO = XgNn'.EQ); 

first excited state 
[El ,mpsl] =minimizeE(hset , D .precision, mpsO) ; 
fprintfC'El = Xg\n',El)j 

As a second example, we focus on the real time evolution with respect to 
the Heisenberg antiferromagnet. Starting state is a product state with all spins 
pointing in z-direction except the central spin which is flipped. We evolve the 
state with the method reduceD and calculate at each step the magnetization 
m z of the central spin. As before, N is equal to 10 and D is set to 5. 

H=10; 
D=5; 
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dt=0.03; 
jf lipped=5; 

'/, magnetization in z-direction 
oset=cell(l,N) ; 

sx=[0, 1; 1,0] ; By=[0,-li;ii,Q]; sz= [1 , ; , -1] ; id=eye{2); 
for j=l:N, oset{l,j}=id; end; 
oset{l, jf lipped}=sz; 

% time evolution operator 
h=kron(sx,sx)+kron(sy,sy)+kron(sz,sz) ; 
w=expm(-li*dt*h) ; 

w=reshape(w, [2,2,2,2]) ; w=permute(w, [1,3,2,4] ) ; w=reshape(w, [4,4]) ; 
[U,S,V]=svd2(w) ; eta=size(S , 1) ; 
U=U*sqrt(S); V=sqrt(S)*V; 

U=reshape(U, [2, 2, eta]) ; U=permute(U, [4,3,2,1] ) ; 
V=reshape(V, [eta, 2, 2]) ; V=permute(V, [1,4,3,2]) ; 
I=reshape(id, [1,1,2,2]) ; 
mpo_even=cell(l,N) ; 
mpo_odd=cell(l,N) ; 

for j = l : N , mpo_even{j }=I ; mpo_odd{j]-=I ; end 

for j=l:2:(N-l), mpo_odd{j}=U; mpo_odd{j+l}=V ; end 

for j=2:2:(N-l), mpo_even{j}=U; mpo_even{j+l}=V; end 

% starting state (one spin flipped) 
mpsO=cell(l,N) ; 
for j=i:H 

if j==jflipped, state=[0; 1]; else state=[l; 0]; end 
mpsO{j}=reshape( state, [1,1,2] ) ; 

% time evolution 
mps=mps0 ; 
mzvalues= [] ; 
for step=l : 50 

fprintf ('Step 7,2d: '.step); 

[mps , K] =reduceD (mps ,mpo_even, D , precision) ; 

[mps ,K] =reduceD (mps ,mpo_odd,D , precision) ; 

mz=expectationvalue (mps , oset) ; 

mzvalues= [mzvalues ,mz] ; 

fprintfCmz^/.gNn'.mz); 

A comparison of the results produced by these examples to exact caclulations 
is shown in figure [26] It can be seen that already for moderate values of D the 
precision is very good. 
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